NASA 

Technical 

Paper 

2221 

January 1984 


AESOP: An Interactive 
Computer Program for the 
Design of Linear Quadratic 
Regulators and Kalman Filters 


i NASA ? L 

| TP 

i 2221 jr- 

i c - 1 ; 



Bruce Lehtinen 
and Lucille C. Geyser 


LOAN COPY: RETURN TO 
AFWL TECHNICAL UBR A T? 
KIRTLAND AFB, MM. 37 il/ 


f\JASA 


TECH LIBRARY KAFB, NM 


NASA 

Technical 

Paper 

2221 

1984 


TECH LIBRARY KAFB, NM 



□ Db7Bt>5 


AESOP: An Interactive 
Computer Program for the 
Design of Linear Quadratic 
Regulators and Kalman Filters 


Bruce Lehtinen 
and Lucille C. Geyser 

Lewis Research Center 
Cleveland, Ohio 


NASA 

National Aeronautics 
and Space Administration 

Scientific and Technical 
Information Office 


1984 



Contents 


Page 

Summary 1 

Introduction 1 

Theoretical Background and Problem Formulation 2 

Open-Loop System Description 2 

Eigenvalues and Eigenvectors 4 

Controllability, Observability, and Mode Shapes 4 

Residues 5 

Steady-State Linear Quadratic Regulator (LQR) Design 5 

Steady-State Kalman Filter Design 7 

Normalization 7 

Stochastic Linear Quadratic Regulator Design 9 

System Response to Noise Inputs 9 

Transfer Functions and Frequency Response Calculations 10 

Transient Response Calculation 13 

AESOP Program Operation 15 

Program Structure 15 

Operating Procedure 15 

Data Input and Output 19 

Description of AESOP Functions 20 

Series 100 - Program Control 22 

Series 200 - Data Input and Revision 22 

Series 300 - Matrix Formation 24 

Series 400 - Open-Loop System Analysis 24 

Series 500 - Frequency Responses and Bode Plots; and Series 700 - Transfer Functions 25 

Series 600 - Transient Responses 26 

Series 800 - LQR and Filter Design 26 

Series 900 - User-Supplied Subroutines 28 

Concluding Remarks 28 

Appendixes 

A -Symbols 30 

B - Subroutine Descriptions 33 

C- Test Cases 64 

Test Case I - Third-Order Problem to Demonstrate Full AESOP Program Capabilities 64 

Terminal Printout for Test Case I 65 

Test Case II - Interactive Design of a Nonzero-Set-Point Regulator 80 

D - Terminal Output Options and Main PROCDEF 102 

Standard Terminal Output 102 

Extended Terminal Output 102 

PROCDEF AESRUN 103 

E - Flow Chart 104 

F - Prerequisite Table 108 

References Ill 


iii 



Summary 

AESOP is a computer program for use in designing 
feedback controls and state estimators for linear multi- 
variable systems. AESOP is meant to be used in an 
interactive manner. Each design task that the program 
performs is assigned a “function” number. The user 
accesses these functions either (1) by inputting a list of 
desired function numbers or (2) by inputting a single 
function number. In the latter case the choice of the 
function will in general depend on the results obtained by 
the previously executed function. 

The most important of the AESOP functions are those 
that design linear quadratic regulators and Kalman 
filters. The user interacts with the program when using 
these design functions by inputting design weighting 
parameters and by viewing graphic displays of designed 
system responses. Supporting functions are provided that 
obtain system transient and frequency responses, transfer 
functions, and covariance matrices. The program can 
also compute open-loop system information such as 
stability (eigenvalues), eigenvectors, controllability, and 
observability. 

The program is written in ANSI-66 Fortran for use on 
an IBM 3033 using TSS 370. Descriptions of all sub- 
routines and results of two test cases are included in the 
appendixes. 


Introduction 

The computer program called AESOP (Algorithms for 
Estimator and OPtimal regulator design) was written to 
solve a number of problems associated with the design of 
controls and state estimators for linear time-invariant 
systems. The systems considered are modeled in state- 
variable form by a set of linear differential and algebraic 
equations with constant coefficients. Two key problems 
solved by AESOP are the linear quadratic regulator 
(LQR) design problem and the steady-state Kalman filter 
design problem. The remainder of AESOP is devoted to 
calculations in support of these two problems, mainly for 
analyzing the open-loop system and evaluating the 
resulting control or estimator designs. Thus the overall 
program can be subdivided as follows: 

(1) Open-loop system analyses 

(2) Control and filter design 

(3) System response calculations 

The AESOP program was developed at Lewis for use 
in conducting design studies in propulsion system control 
(refs. 1 and 2). AESOP was an outgrowth of a previously 
developed control design program called LSOCE (ref. 3), 
which had been used in supersonic inlet controls develop- 
ment (refs. 4 and 5). AESOP differs from LSOCE mainly 
in that it was designed to be operated in an interactive 


manner, whereas LSOCE was strictly a batch type of 
program. In addition, AESOP contains system response 
and evaluation features that are not present in LSOCE. 
These additions tend to enhance AESOP’s use as an 
interactive design tool. 

Other control design computer programs appearing in 
the literature perform computations similar to those of 
AESOP. Notable among the original LQR design 
programs are ASP by Kalman and Englar (ref. 6) and its 
Fortran version VASP (ref. 7). Subsequent LQR design 
packages were the OPTSYS program of Bryson and Hall 
(ref. 8), the ORACLS program of Armstrong (ref. 9), 
and Honeywell’s DIGIKON (ref. 10). Computer-aided 
control system design program development has 
accelerated in recent years. A good summary of this 
development is contained in reference 11. Here, over 20 
control design programs and packages, including 
AESOP, are discussed in varying degrees of detail. They 
represent a variety of design methodologies, ranging 
from classical single-loop approaches to multivariable 
LQR and multivariable frequency domain approaches, 
for both continuous and discrete formulations. Most are 
written in Fortran and have some sort of interactive 
capability, but except for a few commercially available 
packages, most are neither completely documented nor 
generally available. AESOP, at the present time, is the 
only interactive LQR type of control design program that 
is in the public domain. (Contact COSMIC, The 
University of Georgia, Athens, Ga. 30602, concerning 
the availability of this program.) 

The AESOP program is structured around a list of 
predefined and numbered functions. Each function 
performs, basically, a single computation associated with 
control, estimation, or system response determination. 
For example, one AESOP function computes the 
eigenvalues of the open-loop system matrix A, another 
function reads in the A matrix, etc. These functions are 
described fully in the section Description of AESOP 
Functions. The use of these functions and the part they 
play in AESOP can be described in general terms with the 
aid of figure 1 . The figure illustrates what the user of the 
program does (left side of fig. 1), what the program does 
(right side of fig. 1), and the interaction between the user 
and the program. 

The user begins by defining the problem to be solved 
(e.g., by defining the matrices that define the state- 
variable model of the open-loop system). The user then 
provides this information to AESOP as input data, 
generally storing it in a data file. Next the program 
“prompts” the user to enter a list of function numbers 
that are to be performed, in sequence, by AESOP to 
solve the user’s problem. Usually this list of numbers is 
entered at a terminal, but it can also be entered from a 
prestored data file. The AESOP program then executes 
the desired functions in proper sequence, storing away all 
results on an output file (fig. 1) as “off-line output” but 



User 


Program 



Figure 1. -Overview of AESOP program operation. 


also displaying selected portions of the results back at the 
user's terminal (fig. 1) as “on-line display." The user 
then decides whether to terminate the program or to 
request that more functions be performed. The program 
again prompts the user to enter numbers that define the 
new functions, which can be entered singly or as a 
number string. Typically, one of the functions a user 
would enter, at this time, would be one that allows the 
user to vary some problem parameter. In this way the 
user can effectively interact with the program in an on- 
line manner to achieve the desired design results. At the 
conclusion of the terminal session the user commands the 
output data file to be printed and hard copies to be made 
of any graphic output generated that was not previously 
displayed on-line. 

This concludes the overview of the basic operation of 
the AESOP program. The next section describes the 
various design problems that can be solved by using 
AESOP, indicating what function numbers the user 
would request in order to perform the computations. 
Following that section is the section AESOP Program 
Operation. Here, examples are presented of a typical 
dialog between the user and the program. Following that 
is the section Description of AESOP Functions, which 
describes each of the 78 functions, what input each 
requires, and what calculations each performs. These 
latter two sections serve as a guide to which the user can 


refer as necessary while running the program. 
Appendixes include information such as a symbol table, 
brief subroutine descriptions, and two test cases that are 
useful both for program checkout and for gaining an 
understanding of how the program operates. 

Theoretical Background and 
Problem Formulation 

The computations performed by the AESOP program 
can be grouped into five basic categories. These cate- 
gories are illustrated in figure 2. This section presents the 
equations that define the various problems to be solved 
and indicates the solution methods used by AESOP. 
After reading this and the next section, the reader should 
be able to use AESOP to solve a number of “standard" 
problems by using “standard" sequences of function 
numbers. The reader will then be able to devise other 
function number sequences that would allow other more 
specialized problems to be solved. 

Open-Loop System Description 

Before beginning any control or filter design on a 
linear dynamic system, it is important to thoroughly 
analyze the open -loop system under consideration. The 
linear open-loop system defined for use throughout this 
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Figure 2. -Control system design computations performed by AESOP program. 


report is given in the following state-variable form and 
shown schematically in figure 3. The state equation is 

x = A x + B u + D w (1) 

where 

x Nth-order state vector 

u NCth-order control vector 

w NDth-order white-noise disturbance vector 

and A, B, and D are matrices of appropriate dimensions. 
Fortran symbols for matrices used in the AESOP 
program coding are used herein whenever possible. A 
measurement equation, defining the system’s measurable 
output vector is 

z = Z!+v (2a) 

Zj = H x (2b) 


where 

z NMth-order measurement vector 

v NMth-order white-noise measurement vector 

H NM-by-N matrix 

In addition to the measurement vector an output vector, 
which represents unmeasurable outputs, is defined as 

y = C x + DOUT u (3) 

where y is an NOth-order output vector. Finally, a set of 
noise-free measurements called a set-point vector are 
defined. These represent outputs (NC in number) that are 
to be regulated to desired constant set-point values. This 
vector is given as 

y sp = CSPx (4) 

where y sp is an NC th -order set-point vector. 
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Eigenvalues and Eigenvectors 

Of prime importance in designing a control system is 
knowledge of the open-loop system structure and 
stability. This knowledge affects the designer's choice of 
performance index weighting matrices, sensed variables 
to use for control or estimation, etc. 

Open-loop stability is determined by the eigenvalues of 
the system A matrix. The system is stable if and only if 
these eigenvalues X,- (i = 1, 2, . . N) all have negative 
real parts. Consider the unforced version of equation (1), 

x = A x (5) 

Define a new state vector x, relating to x through the 
transformation matrix T as 

T x = x (6) 

Substitute for x in equation (5) to obtain 

x = T~ 1 A T x (7) 

If we let T-i AT be equal to a diagonal matrix A, 
equation (7) can be rewritten as 

x = A x (8) 

The diagonal elements of A are the eigenvalues of A, T is 
the eigenvector matrix (a matrix whose columns are the 
eigenvectors of A), and x is defined as the modal state 
vector. The value of A is computed by using AESOP 
function 501 , and the eigenvectors are obtained by using 
function 402. To avoid complex arithmetic, a block 
diagonal form is used for the matrix A such that a 
complex eigenvalue pair (X/, X/+ j) = (a+y7?* oc-j&) 
appears in the 2-by-2 diagonal block of A, 

a! -0 


0! a 

where the a’s are along the diagonal. Let the complex 
eigenvector pair (77/, 77 / + 1 = y+jb, y-JS) correspond to 
the complex eigenvalue pair a±jP. Then in the columns 
corresponding to the defined diagonal block of A there 
appear two real vectors t, and t/ + 1 defined as 

t/=7 + 5 

t/+l=7-S 


Hence, when A is block diagonal, T is called a modified 
eigenvector matrix. AESOP function 402 also calculates 
the so-called mode shapes . For a real eigenvector the 
mode-shape vector is the same as the eigenvector. 
However, for a complex eigenvector pair the corre- 
sponding mode-shape vector pair contains, in successive 
columns, the magnitude of 77/ and the phase of 777. Each 
mode-shape vector is normalized by dividing all of the 
elements by the magnitude of the largest element. The 
phase of the largest element is set to zero, and the phases 
of all other components of the vector are adjusted 
accordingly. 

Controllability, Observability, and Mode Shapes 

Once the eigenvalues and eigenvectors (mode-shape 
vectors) are calculated, it is an easy task to determine 
controllability and observability . For this purpose, 
system equations (1) and (2a) can be rewritten in terms of 
the modal state vector as 


= Ax + T-l B u 

(9) 

= H T x + v 

(10) 


System controllability is determined by examining the 
elements of T _1 B. The system is uncontrollable if 
elements in a row of T- 1 B are zero, meaning that it is 
impossible to excite a component of the modal state 
vector x with the control vector u. Also, using this modal 
formulation, one can think of the matrix T~ 1 B as being 
the control effectiveness matrix. That is, the relative 
magnitudes of the row elements of T” 1 B define the 
relative influence each control input has on a modal state 
variable (mode). For a meaningful comparison, however, 
the control inputs must be normalized (nondimen- 
sionalized). Normalization can be done by using AESOP 
function number 404. Normalization (and unnormaliza- 
tion) is discussed in detail later in this section. The 
control effectiveness matrix is calculated by AESOP 
function 403. 

System observability can be determined similarly by 
using the modal state form. From equation (10) it can be 
seen that, if all elements of column k of the H T matrix 
are zero, modal state k will be unobservable through 
measurement z. Also, the relative magnitudes of the 
elements of row k of H T define the relative contribution 
each mode makes to measurement z*. This information is 
useful, for example, if one wishes to minimize the 
number of measurements (sensors) required when 
designing a control system that is to shift certain system 
poles (modes, modal states). System observability (the 
H T matrix) is computed in AESOP by using function 
number 403. 
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In AESOP both the controllability matrix (T~1B) and 
the observability matrix (H T) are printed out in mode- 
shape format. This means that, for T-i B, when two 
successive rows k and £+1 relate to a complex modal 
state pair (x*, x* + 1 ), the elements in the columns of 
T~1B are magnitudes and the (A:+l)th elements are 
phase angles. Similarly, for the H T matrix, for a 
complex modal state pair, elements in the column of 
H T are magnitudes and those in the (k+ l)*h column are 
phase angles. 

Residues 

The availability of matrices H T and T-* B makes it 
very easy to compute the system residues. Consider the 
system in modal state vector form given by equations (9) 
and (10). Let B = T-* B and H = H T. Thus equations 
(9) and (10) can be written as 


x = A x + B u 


(10a) 


j 

1 


0 



Substituting from equation (lOf) into equation (lOe), we 
obtain an equation that defines the residue matrices, 
namely. 


N n N U I? n 

i \j y 1 u i d 

s ~h 


OOg) 


z = H x + v 


(10b) 


Thus the y* residue matrix is simply 


For a single-input-single-output linear system a transfer 
function g (s) can be written in so-called residue form as 

N 

*(*>= £ C -V uoo 

j= 1 

where each of the N constants rj is defined as a residue at 
the transfer function pole Xy. The residues define the 
relative magnitude with which the system input affects 
the system output through each system pole. This single 
input/output concept generalizes directly to the multiple 
input /output case. Here the transfer function matrix 
G^) for the system of equations (10a) and (10b) can be 
written as 


Ry=H Ey B (lOh) 

For a real eigenvalue Xy the elements of the corresponding 
residue matrix Ry are real, being computed simply as the 
(outer) product of the y't h column of H and the y th row 
of B. 

For a complex eigenvalue pair (Xy, X y+1 ) AESOP 
makes use of the modified eigenvector matrix form for T, 
which means that H and B are also used in that form. 
Thus real arithmetic can be used in computing the real 
and imaginary parts of the residue matrix. AESOP prints 
out that matrix in polar form (one matrix of residue 
magnitudes followed by one matrix of residue phase 
angles). The residues are computed along with open -loop 
controllability and observability checks in function 403. 


G(5) = H(sI-A) _1 B (lOd) 

or in residue form 

N „ 

G(s)= T. =H(sI-A)- 1 B (lOe) 

H S ~ X J 

where now the N elements Ry are residue matrices . Since 
A is a diagonal matrix, we can rewrite the matrix 
(si — A)- 1 as 

(100 

where 


Steady-State Linear Quadratic 
Regulator (LQR) Design 

One of the primary functions of the AESOP program 
is to compute solutions to the steady-state linear quad- 
ratic regulator problem. Because this problem has been 
well documented (ref. 12, e.g.), the results are only 
briefly summarized herein. The system to be controlled is 
described by 

x = Ax + Bu (11) 

where the state x is assumed to be measurable and no 
plant disturbances are present. 

A control that minimizes the quadratic performance 
index 
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j = f ' “ { X T QC x + 2xT NN u + uT(PCINV) ~ 1 u } dt (12) 
J 0 

is given by 

u = — KC x (13) 

For the optimal solution to exist, weighting matrices QC, 
NN, and PCINV must be as follows: 

(1) PCINV is positive-definite 

(2) QC can be written as QC = MQMT where the 
pair (M,A) is observable and Q is symmetric and 
positive-definite 

(3) QC-NN*PCINV*NN T is nonnegative-definite. 
Feedback gain matrix KC is found by solving the 
following matrix Riccati equation for matrix SS: 

SS(A - B PCINV NNT) + (A - B PCINV NNT) T SS 

- SS(B PCINV BT)SS + (QC - NN PCINV NNT) = 0 

(14) 

Then KC is given as 

KC = PCINV(B T *SS + NNT) (15) 

Figure 4 shows the structure of the LQR solution. The 
gain matrix KC and the Riccati equation solution matrix 
SS are computed in AESOP by function 801. The closed- 
loop state equation for the regulator system shown in 
figure 4 is given by 

x = (A-BKC) x (16) 

AESOP uses the eigenvector decomposition method (ref. 
8) to solve the Riccati equation, and as a byproduct it 
prints out both the eigenvalues and eigenvectors of 

A-B KC. 


The Riccati solution matrix SS theoretically is positive- 
definite and symmetric. Three error checks are provided 
in AESOP (functions 805, 806, and 807) to determine the 
accuracy of the computed SS. The eigenvalues of SS are 
computed and should be positive and real. The 
differences of the off-diagonals are displayed as a 
symmetry check. Finally, the computed SS is substituted 
back into equation (14) and a residual matrix is 
computed. 

The standard steady-state linear quadratic regulator 
problem just outlined assumes that no command inputs 
are present. This problem can be modified to include set- 
point inputs by introducing a set of NC set -point outputs 
defined by 

y sp = CSP x (17) 

These outputs are to be made equal, in steady state, to 
NC corresponding desired set points y sp . This is the 
so-called nonzero-set-point regulator problem of 
Kwakernaak (ref. 12). The solution is to allow a 
feedforward term in the control such that the control law 
of equation (13) is modified to the form 

u = - KC x 4- KFF y spd (18) 

Figure 5 shows the configuration of the nonzero-set- 
point regulator. By stipulating that in steady state 
ysp = yspd » matrix KFF can be computed as 

KFF = [-CSP (A-B KC) 1 B]" 1 (19) 

Thus, with NC degrees of control freedom available, NC 
outputs (y sp ) can be positioned in steady state by using a 
feedforward matrix. The matrix KFF is simply the inverse 
of the closed-loop LQR system transfer function matrix 
evaluated at 5 = 0. 


I 
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Figure 5. - Block diagram of nonzero-set-point linear regulator. 


Steady-State Kalman Filter Design 

The second major computation performed by the 
AESOP program is the design of the steady-state Kalman 
filter for a linear time-invariant system described by 
equations (1) and (2) and shown schematically in figure 3. 
Data required to define the problem consist of plant 
matrices A, B, and H and power spectral density matrices 
for disturbance w and measurement noise v. These white 
zero-mean Gaussian noise signals are described by their 
covariance matrices, namely 

£(w(/)wT (f + r)}=Q<5(r) (20) 

and 

^{vtOvT (t + T)} =(RRINV)“ 1 5(t) (21) 

where matrices Q and (RRINV) 1 are power spectral 
density matrices. For the AESOP program the 
disturbance power spectral density matrix is entered as 
matrix QQ, where QQ is defined as 

QQ = DQDT (22) 

Figure 6 is a block diagram of a linear system in 
“standard” form with its associated Kalman filter. The 
state equation defining the Kalman filter is 

S = (A-KEH) x + Bu + KEz (23) 

The constant gain matrix KE characterizes the filter and 
is obtained in AESOP by using function 809. In ob- 
taining KE, AESOP solves the following Riccati 
equation: 


A PP + PP- AT - PP HT RRINV H PP + QQ = 0 (24) 

Matrix PP is the covariance of estimation error e, where 
e = x-x. The Kalman gain matrix is computed by using 

PP as 

KE = PP HT RMNV (25) 

As is the case with the LQR gain solution, AESOP uses 
the eigenvector decomposition method to solve the 
Riccati equation. Thus as a byproduct the eigenvalues 
and eigenvectors of the Kalman filter (of matrix 
A-KE-H) are printed out. 

As mentioned in the case of the LQR the Riccati 
solution matrix (pp in this case) should be positive- 
definite and symmetric. Three error checks are provided 
in AESOP to check on the accuracy of the estimation 
error covariance matrix PP: 813, to check for positive- 
definiteness; 814, to check symmetry; and 815, to 
perform a residual error check. 


Normalization 

One useful operation that is often performed on the 
matrices appearing in equations (1) to (4), which define 
the open-loop system, is that of normalization. 
Normalization alleviates possible numerical problems; 
allows meaningful comparison between control, state, or 
output variables having different units; and is generally 
recommended for all control and estimator design work. 
Generally one defines a normalization factor (usually a 
full-scale or operating point value) for each component 
of each vector. In AESOP a set of diagonal 
normalization matrices are defined for the system 
variables as follows: 
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The normalization (scale) factors SCX, etc., can be 
considered to be diagonal matrices but are stored in 
AESOP as single-dimensioned arrays. Function 404 is 
provided in AESOP to normalize all of the matrices that 
define the system and the control and estimation 
problems, namely: A, B, C, D, DOUT, CSP, QQ, and 
RRINV. As an example of the calculations performed, 
consider the normalization of the CSP matrix. We have 
that 

y sp = CSP x (26) 

Explicit definition of the normalization factors for x and 
y sp are given by 

y sp £ SC YSP y sp (27) 

x £ SCX x (28) 


CSP = (SCYSP) - 'CSPSCX (30) 

The other matrices are normalized in a similar manner. 
Note that performance index weighting matrices QC, 
NN, and PCINV are not normalized in AESOP because 
they are considered to be “free” parameters to be 
manipulated by the designer. For example, “Bryson’s 
rule,” the often -used rule of thumb for choosing starting 
values of QC and (PCINV)” 1 states that the matrices 
should be diagonal, where each diagonal term is simply 1 
divided by the square of the maximum (or operating 
point) value of the corresponding state or control 
variable. If the system is normalized, the same result can 
be obtained by simply making QC and PCINV identity 
matrices. 

If normalization is used before conducting LQR or 
Kalman filter designs, it may be desirable to have 
normalized gain matrices KC, KE, and KFF put back in 
dimensional form (unnormalized). Function 405 is pro- 
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vided in AESOP for this purpose. In addition, this 
function unnormalizes the error covariance matrix PP. 

Stochastic Linear Quadratic Regulator Design 

The solution of the linear quadratic regulator problem 
requires that the state vector x be completely measurable. 
In general, this will not be possible. Usually, only a 
vector of NM noisy measurements z, which are linearly 
related to the state x, will be present for use by the 
control. In line with the separation principle (ref. 12), the 
optimal control for this situation is constructed by 
feeding back an optimal state estimate (generated by a 
Kalman filter) through the optimal regulator gains KC. 
This system is optimal with respect to minimizing the 
stochastic equivalent of the quadratic performance index 
given by equation (12). That equivalent index is given by 

J=E[ xTQC x + 2xTNN u + uT(PCINV) “ 1 u } (31) 

AESOP provides the means for solving this optimal 
control problem by using the previously mentioned two 
functions for computing gain matrices KC and KE. The 
structure of the complete stochastic LQR problem is 
shown in figure 7. In addition to gain computations it is 
also of interest to compute various system responses to 
characterize the complete closed-loop system. Of 
particular interest in the case of the stochastic LQR 


system are the values of the covariance matrices for the 
system state, control, and output vectors. This is dis- 
cussed in the following section. 

System Response to Noise Inputs 

The primary way to evaluate the overall performance 
of a system controlled by a stochastic linear quadratic 
regulator (such as is shown in fig. 7) is to examine the 
mean square or rms values of the various system vari- 
ables. More generally, the quantities one wishes to 
compute are the covariance matrices for system vectors x, 
x, u, z, and y. In particular, the mean square values are 
the diagonals of the covariance matrices. The two 
covariance matrices are XX, the covariance of the state 
vector x, and PP, the covariance of the Kalman filter 
estimation error e. As was mentioned previously, matrix 
PP is computed by AESOP function 809, which conducts 
the Kalman filter design. The second covariance matrix, 
XX, is obtained by solving the following Lyapunov 
matrix equation (ref. 12): 

(A - B*KC)XX + XX(A - B-KC)T 

+ B KC PP 4 - PP KCT BT + QQ = 0 (32) 

AESOP uses an iterative method developed in refer- 
ences 13 and 14 to solve this equation. In comparison 



Measurement, z 


Figure 7. - Stochastic LQR block diagram showing plant, Kalman filter, LQR gain matrix, and feedforward gain matrix. 
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with other techniques this method has been found to be 
especially effective for cases where system order N is 
large (50 to 100). 

The Lyapunov solution is computed in AESOP 
function 817. Also, this function computes three other 
system covariance matrices, all simple functions of XX 
and PP. They are 

(1) The covariance matrix of control u, 

UU = KC-(XX - PP)*KC T (33) 

(2) The covariance matrix of measurement component 
Zb 

ZZ = H XX HT (34) 

(3) The covariance matrix of output y, 

YY = OPP-CT + (C- DOUTKC) 

(XX - PP) (C - DOUT KC)T (35) 

Note that covariances XX, UU, ZZ, and YY can be 
computed for cases where either (1) no control is used 
(open-loop response) or (2) no Kalman filter is used (state 
feedback only). The open-loop case can be computed by 
simply calling function 817 without first computing either 
PP or KC (these two matrices will thus be all zeros). The 
state feedback case can be computed by first calling 
function 801 to obtain KC and then calling function 817. 

In addition to solving Lyapunov equation (32) for the 
state covariance matrix, AESOP has an error check 
function (number 818) that gives information on the 
accuracy of the solution. Consider a general Lyapunov 
equation 

AX + X AT + W = 0 (36) 

Let the actua^ computed solution to equation (36) be X. 
Substituting X into equation (36) for X we obtain 

AX + X A T + W = R (37) 

where^ R is a residual matrix. Define an error matrix 
E ^X-X. Subtracting equation (36) from equation 
(37), we obtain another Lyapunov equation 

A E + E A T - R = 0 (38) 

AESOP function 818 uses matrix XX obtained from the 
Lyapunov solution of function 817 and (1) solves for the 
residual matrix, (2) solves a Lyapunov equation to obtain 


the error matrix, andJ3) computes an average error value 
as Trace (E)/Trace (X). 

Transfer Functions and Frequency 
Response Calculations 

It is often quite useful to examine the characteristics of 
a state-variable system, either open or closed loop, in the 
frequency domain. For instance, one may wish to analyze 
the pole-zero structure of the system transfer function 
matrix, given a state-variable system description. For this 
purpose, one needs to be able to compute transfer 
function poles, zeros, and gain given the system matrices. 
As another example, one may examine the transfer 
function matrix of an optimal feedback controller to see 
if any simplifying pole-zero cancellations exist such that a 
lower order approximation can be made. Here too it is 
desirable to compute poles and zeros from a state-space 
description. Frequency response plots (for example, Bode 
plots) may be desirable so that one can evaluate, using 
classical frequency domain criteria, the response of a 
control system that was designed by using LQR methods. 
Or, given a state-space, open-loop system description, 
one may wish to compute a matrix of system transfer 
functions or a matrix of frequency responses that can 
subsequently be used by a frequency domain control 
design program. For these reasons, it was decided to 
include in AESOP the capability to compute transfer 
functions and frequency responses for various systems 
and subsystems defined in state-space terms. 

Consider the generalized rtt h -order system described by 
state equations 

x = Ax + Bu (39) 

y = Cx + Du (40) 

and having “nc” inputs u and “no” outputs y. These 
equations could represent any linear system (open-loop 
plant, Kalman filter, closed-loop regulator, etc.) by 
appropriate choice of vectors x, y, and u and matrices A, 
B, C, and D. The transfer function matrix G(s) relating 
output vector y to input vector u can be written as 

y(s) =[C(sI- A) -1 B + D]u(s) =G(s)u(s) (41) 

AESOP allows the user to obtain solutions to equation 
(41) for a variety of system configurations. It computes a 
transfer function G y(s) relating a component u j(s) to a 
component y ,*(s) in two forms. 

The first transfer function form computed is where 
each Gy (s) is a ratio of polynomials. In this case the 
general expression for Gy(s) is 
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Gy (5) = 


n + 1 

£ fl / t - i **- 1 

Ar= 1 


5«+ £ b k - 1 S*- 1 
Ar= 1 
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Note that since a feedforward matrix D is assumed, 
Gij(s) may have equal -order ( n ) numerator and 
denominator. AESOP uses the technique of reference 15, 
modified slightly to include a D matrix, to compute 
coefficients a % and b Constants b k are the coefficients 
of the characteristic equation of A and are thus common 
to all G ij(s). 

The second transfer function form computed by 
AESOP is the so-called factored form, where the general 
expression is given as 

m 

Kijjl(s+z k ) 

Gij(s) = —f (43) 

^ S+Pk) 


where Zk and pk are the transfer function zeros and poles, 
respectively. In addition to poles and zeros, AESOP 
computes gains Ky and the number of numerator zeros 
m (m <ri). The poles Pk are obtained in AESOP simply as 
the eigenvalues of A. The method for obtaining m and z k 
depends on the value of D,y. 

For cases where D z y is equal to zero the values of m and 
Zk are obtained by using a method developed by Davison 
(ref. 16). The method is essentially based on a concept 
from root locus theory. That is, if a proportional loop is 
closed between output y / and input fly and the loop gain is 
allowed to increase to infinity, m of the root loci (poles of 
the closed-loop transfer function) will go to the m open- 
loop transfer function zeros, and the remainder (n~m) 
will go off to infinity. Davison’s method successively 
computes the eigenvalues of such a system while 
increasing the loop gain. It stops when the n — m 
“extraneous” eigenvalues all exceed a “large” value. 

For cases where the value of D z y is nonzero, the number 
of zeroes and poles of the transfer function Gy(s) are 
both equal to n. In this case the zeroes are simply the 
eigenvalues of the matrix 



where 

B ^ [b, ! b 2 : : 


c 


T 

no 


This fact can be seen by applying a feedback control 


Uy= -Ar*y/ 


(45) 


to the system of equations (39) and (40) and allowing gain 
k* to go to infinity. In the limit the eigenvalues of A* 
become the zeros of the transfer function Gy (s). 

The remaining transfer function term in equation (43) 
to be computed is gain Ky . AESOP uses the technique 
described by Brockett (ref. 17) to compute this gain as 


Kij=< 


d //> 


C.*A n-m- \ by, d z y = 0 




(46) 


Computations for two types of transfer functions, 
polynomial form and factored form, are performed in 
AESOP series 500 and 700, respectively. System 
configurations for which calculations are done are (1) an 
open-loop system, (2) a system with state feedback, (3) a 
system with a Kalman filter in the feedback loop, and 
(4) an optimal controller. Polynomial-form transfer 
function coefficients are computed for the purpose of 
plotting frequency responses. The AESOP program uses 
the same subroutines to compute frequency responses 
and transfer functions for each configuration, first 
forming the appropriate A, B, C, and D matrices as 
functions of A, B, KE, KC, etc. Factored-form transfer 
function information (poles, zeros, and gain) is mainly of 
interest in one of two instances: (1) when investigating the 
structure of the open-loop system, and (2) when 
examining the pole-zero structure of an optimal 
controller. Therefore data are obtained by AESOP only 
for these two configurations. Frequency responses, 
however, are obtained for all four configurations 
mentioned previously. Table VII in the next section 
outlines in detail which calculations AESOP does 
perform. 

Figure 8 shows the open-loop configuration for which 
transfer functions and frequency responses can be 
calculated. Corresponding to the form of equation (41), 
the four transfer functions that AESOP computes here 
are 
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Figure 8. - Block diagram of open-loop system for which transfer functions and frequency responses are obtained. Inputs are u and w; outputs 
are Z\ and y. 


Control to output: 

y (5) = t(5l - A) - 1 B + DOUT]u (s) (47) 

Disturbance to output: 

y(s) =C(sI- A) -1 D w(s) (48) 

Control to measurement: 

z (s) = H (si - A) " 1 Bu(s) (49) 

Disturbance to measurement: 

z(s)=H(sI-A)" 1 D w(s) (50) 


Figure 9 shows the second configuration -a system 
with state-variable feedback. Here AESOP computes the 
following frequency responses by using disturbance w as 
the input: 

Disturbance to output: 


y(s) = (C-DOUT-KC)[(sI — (A-B KC)] _1 D w(s) 

(51) 

Disturbance to control: 

u (s) = - KC[sI - (A - B-KC)] “ 1 D w (s) (52) 

Disturbance w is not explicitly used in the design of the 
(noise free) linear quadratic regulator (e.g., eq. (14)). 
However, it is instructive to examine its disturbance 
response in the frequency domain. 

Figure 10 depicts the configuration where disturbance 
w is explicitly considered in the control system design, 
namely the stochastic linear quadratic regulator. The 
overall system order in this configuration is 2N, there 
being N states x and N state estimates x. In obtaining 
frequency responses AESOP first forms partitioned sys- 
tem matrices before calling the subroutine that computes 
the responses. The following responses are computed: 

Disturbance to output: 


Control, u 



Figure 9. -Block diagram of LQR configuration used for obtaining frequency responses input is w; outputs are>> and u . 
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I Kalman filter j 

Figure 10. - Block diagram of stochastic LQR configuration used for obtaining frequency responses. Input is w; outputs are y , z, and u. 


y(s)= CTOT (si - ATOT) ~ 1 DTOT w (s) (53) 

Disturbance to measurement: 

z(s) = HTOT(sI — ATOT) -1 DTOT w(s) (54) 

Disturbance to control: 

u(s) = KCTOT (si - ATOT) ' 1 DTOT w(s) (55) 

where 

[ A ! -BKC "I 
1 , 2N x 2N 

KE H | A-B KC-KE HJ 

r°i 

DTOT = I , 2NXND 

CTOT = [C | -DOUTKC], NOx2N 
HTOT = [H 10], NMx2N 


Note that measurement noise, although a consideration 
in the Kalman filter design, is not considered here as an 
input. 

The last configuration for which frequency responses 
as well as transfer function zeros and gain is computed is 
the optimal controller, shown in figure 11. This is simply 
the feedback portion of the control loop of figure 10, the 
stochastic regulator. The desired response is for measure- 
ment z as an input and control u as an output. This 
transfer function can be expressed as 

u(s) = — KC[sI- (A — B-KC — KEH)] - 1 KE z(s) (56) 

It has been noted in the literature that this optimal 
controller transfer function can have some interesting 
properties; for example it can sometimes be unstable (ref. 
18) or can have right-half-plane zeros (ref. 1). 

Transient Response Calculation 

The AESOP program computes and plots transient 
responses for several important system configurations. 
Consider again the general time-invariant system, which 
is described in state-variable form as 

x(/)=Ax(0 +B u(0 (57) 


KCTOT = [0 | -KC], NCx2N y (0 =C x(t) +D u(/) (58) 
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Kalman filter 


J 


Figure 11. - Optimal controller, comprising Kalman filter and LQR gain matrix. Input is measurement Z and output is control u. 


In AESOP it was desired to solve these general equations 
(57) and (58) for two cases: (1) an initial condition x(0), 
and (2) a step change in input u(/). The normal 
approach, and that taken in AESOP, is to discretize 
equation (57) so as to obtain the exact solution at time 
points to, /i, . . tk that are spaced DT seconds apart. 
The resulting discrete difference equations are thus 


*n+l (59) 

y n = Cx n + f)u n (60) 

where 

*n + l=*(/„ + i) (61) 

Xo = x(/ 0 ) = x(0) (62) 


and matrices # and T are given by the series 
representations 


(1) Open-loop system — The equations for this system, 
shown in figure 3, are 

x=Ax+Bu (65) 

y = C x + DOUT u (66) 

AESOP calculates and plots x and y for ITRMX time 
points for either 

(a) A step input applied to component i of input 

u(0 

(b) An initial condition x(0) on the /th component of 
the state vector x(Q. 

These computations are performed by AESOP functions 
601 and 602, respectively. 

(2) Closed-loop linear quadratic regulator - Figure 4 
shows the configuration of this system. The describing 
equations are 

x = (A — B'KC) x (67) 


4>(DT) = e(ADT) = I + A-DT + 


A2-DT2 

2 ! 


+ . . . 


(63) 


and 


■■ [CH* 


=(■• 


^ A-DT2 A2-DT3 

DT + — — — + — — — + . 


2! 


3! 


•) 


B (64) 


y = Cx + DOUT u (68) 

Matrices A, B, C, and D are appropriately formed by 
AESOP and variables x, y, and u = - KC x are calculated 
and plotted. This is done for an initial condition X/(0) on 
a selected component of x(0). Function 603 accomplishes 
these calculations. The responses for x, y, and u can be 
compared by the user with the open-loop responses in 
order to help evaluate a particular LQR design. 

(3) Nonzero-set-point linear quadratic regulator - This 
system, shown in figure 5, is described by one state and 
three “output” equations: 


Matrix T is computed by assuming that input u (t) is 
constant over t n <t<t n + \. The series in equations (63) 
and (64) are carried out until the desired accuracy in $ 
and T is achieved. 

The general procedure just described is used to 
compute transients for the following three situations: 


x = (A-B KC) x + B-KFF y spd (69) 

y = (C - DOUT KC) x + DOUT KFF y spd (70) 

y sp = CSPx (71) 

u = — KC x + KFF y spd (72) 
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Here again, AESOP forms the appropriate matrices 
before computing the responses. The response of interest 
here is to a step change in a selected component of the set- 
point vector y sp d- Of interest to the user in this con- 
figuration are such things as rise time and overshoot of 
the component of y sp in response to a step change in 
the fcth component of y spc j; interaction - the response of 
the £ th component of y sp to a change in the / th component 
of y spc j; and the magnitudes of control variable 
excursions. The computation and response plotting for 
the nonzero-set-point regulator are performed by 
AESOP function 604. 

AESOP Program Operation 

This section provides an overview of how to operate 
the AESOP program. In particular, it discusses how the 
user interfaces with the program within the IBM 370 TSS 
PCS (Program Control System) environment, how to 
enter data, how to enter a string of control numbers that 
dictate the series of AESOP functions to be performed, 
and in general, how the user interfaces with the program 
at the “terminal,’ ’ not Fortran, level. (It is assumed that 
the reader is familiar with PCS commands and the 
operation of the 370 TSS system.) 

Program Structure 

The general structure and operation of AESOP was 
described in the Introduction and shown schematically in 
figure 1. A more detailed diagram of the program show- 
ing the main AESOP program and nine main subroutines 
is given in figure 12. Each main subroutine contains a 
number of related AESOP functions. The numbers above 
the main subroutine boxes in figure 12 indicate the 
function series that each main subroutine contains (800 
contains functions 801 to 899, etc.). The details as to 
what each function does will be discussed later in this 
section. 


All input and output performed by the main AESOP 
program relate to program control and required inter- 
facing with the user. The main program (1) initializes 
default or reference values, (2) accepts a string of 
function numbers that the user wishes to have performed, 
(3) performs checks to see whether the user has requested 
functions to be done in a reasonable order, and (4) calls 
the appropriate main subroutines in which the desired 
functions reside. 

To run the AESOP program, the user First calls the 
PROCDEF (for PROcedure DEFinition) AESRUN, 
defined in appendix D. Through use of TSS/370 datadef 
(DDEF) statements, AESRUN 

(1) Defines (“datadefs”) the library dataset (file) that 
contains the compiled AESOP program and all sub- 
routines and then loads the program 

(2) Defines the link to the graphics package that 
contains graphics subroutines called by AESOP 

(3) Links (“datadefs”) Fortran unit numbers with 
specific dataset (file) names so that these files can be 
written to or read from during the subsequent AESOP 
run 

AESRUN has a single parameter that is used to identify 
all output datasets to be generated during the AESOP 
run. After calling AESRUN, the user types “AESOP” to 
run the program. 

Operating Procedure 

The key element in the use of AESOP is the single- 
dimensioned function number array (Fortran symbol 
IFN). The user loads, in sequence, this vector with the 
numbers of the AESOP functions that are to be per- 
formed. The use of the IFN vector and the general 
operation of AESOP will be explained with the aid of the 
flow chart of figure 13. The user many also refer to the 
terminal listing included in appendix C for test case I for 
an actual example of how AESOP is run. First, the user 
types “AESRUN” followed by its parameter. The 



Figure 12. - AESOP program structure. 
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AESRUN [I. D. PARAMETER] 


| CALL t 

*ES0P 




I IONPL 
|lPRT = 

T = 0 

1 



Heavy blocks indicate 
occurrence of prompts or 
terminal responses 


r®C 


DO YOU WISH TO MAKE 
PLOTS, Y OR N? 


')©■ 


ENTER THE PLOT NAME - 
8 ALPHANUMERIC CHARACTERS! 


r®C 


DO YOU WISH TO MAKE 
ON-LINE PLOTS. YOR N? 




IONPLT ■ 1 


r©C EXTENDED TERMINAL OUTPUT? ^Yes) — IPRT - 2 — | 

pg) ( READ IN N1 FROM STORAGE? ] Q-* | DATASCT^h 


■ENTER TODAY'S DATE (LESS THAN 
|OR EQUAL TO 20 CHARACTERS) J 

J 


1 ENTER THE PARAMETER YOU 
| USED FOR AESRUN 




f ♦ 


EXTENDED TERMINAL OUTPUT? 


XSH 


IPRT - 2 


TO THE 
DATASET 


4 ENTER NAMELIST N1 DATA (IFN) 

1 i : 

MZ = 0 


I 


Figure 13. -Flow chart for AESOP program. 


computer responds with information that the libraries 
have been set up and all datasets have been datadeffed. 
The user then types “ AESOP’ ’ to call the main program. 
(Note that the heavy blocks in figure 13 indicate either 
user input commands or messages printed out at the 
user’s terminal.) The program then responds by prompt- 
ing the user with the message: 

DO YOU WISH TO MAKE PLOTS, Y OR N? 

If the response is “Y,” the user is then asked to 

ENTER THE PLOT NAME - 8 ALPHANUMERIC 
CHARACTERS 


This name will be used to name the plot dataset in which 
any plots generated will be stored. Next, the user is asked, 

DO YOU WISH TO MAKE ON-LINE PLOTS, Y or N? 

If the reply is “Y,” the program sets an internal flag that 
will cause the program to PAUSE after each on-line plot 
is displayed to allow the user to view it. Off-line plotting 
requires no such action to be taken. Finally, the user is 
asked for two pieces of information that will appear on 
all plots for identification purposes: 

ENTER TODAY’S DATE (LESS THAN OR EQUAL 
TO 20 CHARACTERS) 
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1 


MZ = MZ + U 


Execute function 
IFN(MZ) 


<DC IFN(MZ):999 D©~] 


-0C IFN(MZ):100 


oc IFN(MZ);0 DG> 


~(no) ^Is IFN(JVIZ) the number of an executable FCN?^ Yes) — | 
Prerequisite 


error 


Prerequisite 

check 


No error 


YOU HAVE NOT EXEC- 
UTED THE FOLLOWING 
PREREQUISITE 
jFUNCTIONS 


THESE ARE THE FUNCTIONS YOU HAVE RUN THUS 
FAR. THE PRESENT FUNCTION IS NOT A LEGITI- 
MATE ONE. YOU WILL NOW HAVE A CHANCE TO 
REPLACE IT WITH ONE GOOD ONE, TO TERMINATE, 
OR TO CENTER A LI ST OF FUNCTIONS. 

TYPE Y TO ENTER ONE REPLACEMENT FUNCTION; 
TYPE N TO ENTER A LI ST OF FUNCTIONS 


I @ C Reply = ? ) (Ye§) 1 


TO COMPUTE FURTHER, 
ENTER IFN03), ONE PER 
LINE; TO TERMINATE, 
ENTER 999 (LAST ENTRY 
MUST BE A RETURN) 


■TYPE IN AN 13 NUMBER 
■TO REPLACE THE CURRENT 
JUNCTION; IF THE NUMBER) 
IS 999, THE PROGRAM 
■WILL TERMINATE 


Read new function 
number into IFN 


Read IFN(MZ) 


^ 

I — ©C STORE THE N1 FOR THIS RUN? 


DDEF NAME OF N1 
DATASET to 21 


2-1 


End 


Figure 13. -Concluded. 


and 

ENTER THE PARAMETER YOU USED FOR 
AESRUN 

In view of the fact that graphics subroutines are rather 
specific to the user’s computer system, they are not 


provided as part of AESOP. Thus the preceding messages 
and prompts would be tailored by each user to reflect the 
specific graphics package being used. 

After the graphics-oriented prompts are handled, the 
program displays the message 


EXTENDED TERMINAL OUTPUT? 
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(This question and all subsequent ones are to be answered 
“Y” (yes) or “N” (no). All output produced by the 
program is stored on the output dataset, which the 
PROCDEF AESRUN datadeffed to unit 06. Two subsets 
of this output are available for display on-line at the 
user’s terminal. The EXTENDED TERMINAL 
OUTPUT option (appendix D) allows the user to view a 
more extensive subset of data if desired. 

Next, the user is prompted with the message 

READ IN N1 FROM STORAGE? 

The N1 referred to here is the name of a Fortran 
NAMELIST that contains the vector IFN. For problems 
that are solved over and over again, using the same 
AESOP functions but different data, it is convenient to 
store the NAMELIST N1 in a dataset to avoid having to 
type it each time at the terminal. If this is the case, the 
user would respond with “Y” and the program would 
print out 

DDEF 21 TO THE N1 DATASET 

At this point the user would then type the required 
datadef statement, for example, 

DDEF FT21F001,VS,MYN1DS 

where dataset MYN1DS would contain typical 
NAMELIST data such as 

&N1 IFN = 201,701,999 &END 

However, to enter a new set of numbers into the IFN 
vector, the user would enter “N” and the program would 
respond with 

ENTER NAMELIST DATA AS ’&N1 IFN = , , , &END 

At this point the user would type a NAMELIST Nl, such 
as was given in the preceding example. 

The IFN string can be terminated in two ways. If the 
user ends it with a number greater than or equal to 999 (as 
was done in the preceding example), the AESOP program 
will execute all of the requested functions and then 
terminate the run. (The number 999 is a request for 
termination.) If the user ends the string with the number 
of an executable function, the program will allow more 
function numbers to be entered when the present string of 
functions has been executed. In either case, as soon as the 
NAMELIST has been read by the program, the program 
displays all of the elements of the IFN vector at the 
terminal and starts to execute the requested series of 
functions. 


Figure 13 shows that the main program indexes integer 
MZ each time just before it begins to execute a function. 
This integer denotes the element of the IFN vector that 
contains the number of the function about to be 
executed. The program then performs four successive 
checks on IFN (MZ), the MZd* element of IFN. These 
checks are 

(1) Is the IFN(MZ)>999? If so, terminate the 
program; if not, continue checking. 

(2) Is IFN(MZ) < 100 but not equal to zero? If so, you 
have requested a nonexisting function and will be allowed 
to change the requested number. 

(3) Is IFN(MZ) = 0? If so, you have reached the end of 
the function number string requested and will now be 
able to enter additional function numbers if desired. 

(4) Is IFN(MZ) the number of an existing executable 
function in series 100 to 900? If so, the function will be 
executed; if not, you will be allowed (as in check (2)) to 
change the requested number. 

As an example, first, assume that the present IFN(MZ) 
encountered is a nonexisting function number (i.e., 
checks (2) or (4) above were failed). The program dis- 
plays the messages 

THESE ARE THE FUNCTIONS YOU HAVE RUN 

THUS FAR 

nnn 

THE PRESENT FUNCTION IS NOT A LEGITIMATE 
ONE, YOU WILL NOW HAVE A CHANCE TO 
REPLACE IT WITH ONE GOOD ONE OR TO 
TERMINATE OR TO ENTER A LIST OF FUNC- 
TIONS. TYPE Y TO ENTER ONE REPLACEMENT 
FUNCTION; TYPE N TO ENTER A LIST OF 
FUNCTIONS 

If the user types “Y,” the program responds with 

TYPE IN AN 13 NUMBER TO REPLACE THE 
CURRENT FUNCTION; IF THE NUMBER IS 999, 
THE PROGRAM WILL TERMINATE 

The user would then enter a new function number and the 
program would repeat the four checks. If the new 
number is that of an executable function, the program 
proceeds to execute the function. 

If the user types “N” in response to the previous 
prompt, the program will display the message 

TO COMPUTE FURTHER, ENTER NEXT 
FUNCTION NOS. (13), ONE PER LINE; 

TO TERMINATE ENTER 999 (LAST ENTRY MUST 
BE A RETURN). 

The user may then enter a series of function numbers, 
hitting “RETURN” after each three-digit number and 
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hitting ‘‘RETURN” twice after entering the last number 
in the series. The program also responds with the 
preceding prompt whenever it encounters IFN(MZ) = 0, 
which is the indication that the end of the previously 
requested function string has been reached. 

The previous paragraphs have dealt with situations 
where the program detects nonexistent function numbers 
in the input string. In the case where the function number 
is correct, the main AESOP program then calls the 
appropriate main AESOP subroutine in which the func- 
tion resides. At that point a prerequisite check is begun. 
As a user aid the AESOP program contains a table of 
prerequisite functions (appendix F). The program checks 
the prerequisite table immediately before executing each 
requested function to insure that the specified pre- 
requisite functions have already been performed. For 
example, suppose the user enters function number 401, 
which requests that open-loop eigenvalues (of the A 
matrix) be computed, without preceding this number 
with either number 201 or 202, the functions that form or 
read in data that define the system. Just before calling 
function 401 the program checks the prerequisite table 
and detects that the proper prerequisites (either function 
201 or 202) have not been performed prior to this func- 
tion. The program then displays the message 

YOU HAVE NOT EXECUTED THE FOLLOWING 
PREREQUISITE FUNCTION(S) FOR FUNCTION 401 

and then prints out pertinent function numbers. It then 
displays the message 

IF YOU THINK YOU KNOW WHAT YOU ARE 
DOING AND WISH TO IGNORE THE PREREQS 
AND CONTINUE ON TO DO THIS FUNCTION, 
TYPE Y AND RETURN; OTHERWISE, JUST 
RETURN 

The prerequisites in the table for each AESOP function 
were selected so as to catch most, if not all, errors a user 
might make in selecting a series of interdependent func- 
tions that could lead to calculations that would produce 
major errors. These would be errors such as zero divides 
during inversion of a singular matrix, etc. It was felt that 
protecting against these nuisance errors would make it 
less likely that a user would have to restart the program in 
midcourse because of a major nonrecoverable error. 
However, it is still possible to select a series of functions 
that produce nonsensical results even though prerequisite 
checks show no error. 

To terminate the program, the user types “999” when 
the program asks for more function numbers. The pro- 
gram will then display the message 

STORE THE N1 FOR THIS RUN? 


This allows the user to store the IFN vector, which was 
just used in the present run, on a dataset for possible 
future use. If this is desired, the user replies with “Y” 
and will receive the message 

DDEF NAME OF N1 DATASET TO 22 

The user would execute the required datadef, using 
whatever dataset name was desired, and then type “GO” 
to cause the program to terminate. If the IFN array is not 
to be stored, a previous response of “N” would 
terminate the run immediately. 

Data Input and Output 

A key aspect of the AESOP program is the handling of 
data input and output. The primary input and output 
device is the user’s terminal. Most of the data sent or 
received through the terminal pertain to program control, 
as was demonstrated in the examples in the previous 
section. That section also showed two examples of how 
the program accesses data that are stored in a dataset 
(IFN, read from unit 21) and how the program writes 
data out to a dataset (writing IFN onto unit 22). In the 
program, however, 14 unit numbers are used for input 
and output of data (in addition to unit 02, reserved for 
the terminal, and unit 06, reserved for the high-speed 
printer output). All units are listed in table I together with 
the names, contents, and format of their associated 
datasets. Dataset names are created by the PROCDEF 
AESRUN, which appends the characters in parameter $1 
to an identifying prefix. For example, referring to table I, 
if the parameter $1 were entered as XYZ, the PROCDEF 
would datadef unit 08 to dataset CGXYZ. This dataset is 
the one in which the control gain matrix KC is stored. 
Other datasets that are identified by using the parameter 
$1 are those datadef fed to unit 09 (containing the Kalman 
filter gain matrix), units 10 to 14 (containing frequency 
response magnitudes and phase angles), unit 15 
(estimation error covariance matrix), unit 16 (control 
Riccati equation solution matrix), and unit 17 
(feedforward gain matrix). The remaining units, 33 and 
34, are for datasets that store the data that define the 
design problem to be solved by AESOP, namely the 
open-loop plant data (unit 33) and normalizing factors 
(unit 34). The user must datadef these units and specify 
the names of their associated datasets. 

Much of the input data for AESOP is in NAMELIST 
form. Table II lists all NAMELIST’S used by the 
program and the names of the Fortran variables 
contained in each. A key NAMELIST is MATDAT, 
which contains all matrix coefficients and dimensions 
needed to define the problem to be solved. NAMELIST’S 
CONPAR and ESTPAR, which are useful when 
modifying basic problem data, contain subsets of the 
variables contained in MATDAT. CONPAR is used for 
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TABLE I. - DATASETS AND UNITS USED FOR PROGRAM INPUT AND OUTPUT 


Unit 

Dataset 

Contents of dataset 

f Format 

Function where 

Comments 


name 




read /written 


02 


All terminal Input and output 

Various 

Many 

Terminal input and output 

06 

0UT?l a 

High-speed-pr inter output 

Various 

Many 


08 

CGS1 

KC, control gain matrix 

Unformatted 

205/802 


09 

EG£1 

KE, Kalman filter gain matrix 



206/810 


10 

PFRUZ21 

Frequency response, z( IMEAS) /u( JINC) 




See table VI for details 







on frequency responses 

11 

PFRUY21 

Frequency response, y( I0UT)/u( JINC) 





12 

PFRWZ31 

Frequency response, z( IMEAS) /w(JIND) 





13 

PFRWYS1 

Frequency response, y( I0UT)/w( JIND) 





14 

CFR£1 

Frequency response, u( JINC) /z( IMEAS) 





15 

PPS1 

PP, estimation error covariance matrix 



208/816 


16 

ssn 

SS, control Riccati solution matrix 



207/808 


17 

FFG21 

KFF, feedforward gain matrix 



209/820 


21 

(b) 

IFN vector, contained in NAMELIST N1 

NAMELIST 

Main program 

For input of IFN 

22 

(b) 

IFN vector, contained in NAMELIST N1 



Main program 

For output of IFN 

33 

(b) 

Open-loop plant data, NAMELIST' s MATDAT and REFS 



202 

Usual source of problem data 

34 

(b) 

Normalizing factors; NAMELIST NRMS 



404 and 405 



a $l is the parameter for PROCDEF AESRUN; it serves as an arbitrary identifying tag (fewer than four characters) for the 
naming of datasets), 
b User specified. 


TABLE II. - NAMELISTS USED IN AESOP PROGRAM 


NAMELIST 

Variables in NAMELIST 

Input source 

Comments 

CONPAR 

QC, NN, PCINV 

Terminal 

For revising weighting matrices; 
read in function 203 

ESTPAR 

QQ, RRINV 

Terminal 

For revising noise power spectral 
densities; read in function 204 

MATDAT 

A, 8, C, 0, H, D0UT, CSP, 
QC, NN, PCINV, QQ, RRINV 
N, NM, NC, ND, NO 

Terminal or dataset 

Primary means of entering system 
data; read in function 202; 
changed in function 210 

NRMS 

SCX, SCU, SCY, SCZ, SCYSP 

Dataset 

Contains normalizing factors; 
read in function 404 or 405 

N1 

IFN 

Terminal 

Program control; read or written 
in main program 

REFS 

TSFTR, DT, FI, DELF, ZERMAX, 
AMPSP, AMPSR, AMPICX, IF, ISPACE, 
I OUT, IMEAS, JINC, JIND, ITRMX, 
NCURV, LINL0G, MSPY, MSPYSP, 

MSPU, MSROLY, MSROLX, MICCLY, 

Terminal or dataset 

See table III for definition of 
all variables and default values; 
set in main program; change in 
function 101; read in function 
202 


MICCLX, MICCLU, MICOLY, MICOLX 




making changes in performance index weights, and 
ESTPAR is used for varying noise matrix elements. 

NAMELIST NRMS contains normalizing factors 
whose use is described by equations (26) to (30). The 
remaining NAMELIST, REFS, is used for setting various 
parameter values that specify such things as time steps, 
frequency point spacing and numbers, perturbation 
amplitudes, and input and output selection indices. Table 
III specifies each parameter in REFS and its default value 
and indicates in which AESOP function each parameter 
is used. For a more detailed explanation of each 
parameter, refer to the appropriate function description 
given in the next section. 


Description of AESOP Functions 

Each AESOP function will now be described in 
sufficient detail so that this section can serve as a primary 
reference when using the program. All functions, as 
indicated previously, are grouped into nine series (series 
100 to 900). When possible, a general description will be 
provided for each series, with an accompanying table 
included to show similarities and differences among 
functions within the series. Also, whenever possible, the 
reader will be referred to appendix C, test case I, for an 
example of the use of each function. 
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TABLE III. - DEFINITION OF PARAMETERS CONTAINED IN NAMELIST REFS 


Variable 

Dimension 

Definition 

Default value 

Functions where used 

DT 

Scalar 

Time step 

0.05 sec 

Series 600 (transient responses); 





functions 601 to 604 

ITRMX 

Scalar 

Desired maximum number of time steps 

100 



AttPSR 

NCMAX 

Open-loop system step input amplitude vector 

All elements are set to 1.0 



AMPSP 

NCMAX 

Closed-loop system step input amplitude vector 

All elements are set to 1.0 



AMPICX 

NMAX 

Initial-condition amplitude vector 

All elements are set to 1.0 



MSR0LY a 

NCMAX b , NOMAXC 

Input/output selection matrix 

All elements are set to 1 



MSROLX 3 

NCMAX b , NMAX C 

i 

1 




MICOLY 3 NMAX b , NOMAX c 

MICOLX® NMAX b , NMAX C 

M ICCLY® NMAX b , NOMAX c 

MICCLX® NMAX b , NMAX C 

MICCLU 3 NMAX b NCMAX C 

MSPYSpa NCMAX b , NCMAX C 


MSPY a 

NCMAX b , NOMAXC 





MSPU a 

NCMAX b , NCMAXC 

i 

V 


f 

FI 

Scalar 

Starting frequency for frequency responses 

0.01 Hz 

Series 

500, frequency responses 

DELF 



Frequency-point spacing 

0.02 Hz 




NCURV 



= 2, cross plot open- and closed-loop responses; 

2 







= 1, closed loop only 





L INLOG 



= 1, linear Bode plots; = 2, log Bode plots; 

3 







* 3, both linear and log plots 





IF 



Desired number of frequency response points 

49 




ISPACE 



Every ISPACE frequency response point is listed 

1 




TSFTR 



Time scale factor, used for increased precision 

1.0 




IOUT 

Scalar 

Index for selecting component of y vector 

1 

Series 

500 (frequency responses) 






and series 700 (transfer functions) 

IMEAS 



Index for selecting component of z vector 

1 




JINC 



Index for selecting component of u vector 

1 




JIND 



Index for selecting component of w vector 

1 




ZERMAX 

Scalar 

10 times value of largest expected transfer 

100 rad/sec 

Series 

700 




function zero 





a See table VII for detailed description. 
b Row size. 
c Column size. 



■Ill I III III 


When running the AESOP program, it has been found 
helpful to have available a brief list of all possible 
functions. This list is provided in table IV. Once the user 
has read this report and is generally familiar with what 
the program can do, table IV can be used as a ready 
reference guide while running AESOP. More specific 
details on the various functions are given here. 

Series 100 - Program Control 

Series 100 contains only two functions: one for 
changing reference parameter values, and the other to 
allow the user to change datadefs in the middle of a run. 

101 - Change reference values by using NAMELIST 
REFS. -Table III describes various parameters in 
NAMELIST REFS that are used in the AESOP program 
to determine time and frequency steps, input and output 
indices, etc. The table also gives the default values of 
these parameter values that are initialized in the main 
AESOP program. Function 101 allows the user to 
change any or all of these parameters. It prompts the user 
with the message 

ENTER CHANGES TO NAMELIST REFS (TSFTR, 
DT, FI, DELF, ZERMAX, AMPSP, AMPSR, 
AMPICX, IF, ISP ACE, IOUT, IMEAS, JINC, JIND, 
ITRMX, NCURV, LINLOG, MSPY, MSPYSP, MSPU, 
MSROLY, MSROLX, MICCLY, MICCLX, MICCLU, 
MICOLY, MICOLX) 

after which the user can enter the desired parameter 
changes. An example of the use of function 101 appears 
in appendix C, test case I, page 65. 

102— PA USE to allow user to change datadefs. - If the 
user requests this function, the program will display the 
message 

YOU MAY NOW CHANGE YOUR DDEFS IF YOU 
WISH. DON’T FORGET TO CLOSE AND RELEASE 
THE OLD ONES FIRST 

The program then effects a Fortran PAUSE and the 
keyboard unlocks to allow the user to make the appro- 
priate changes. This function is useful, for example, 
when the user wishes to read in NAMELIST's MATDAT 
and REFS from dataset BBB, where previously similar 
data had been read from dataset AAA. At the requested 
PAUSE, the user can then enter 

CLOSE AAA 
RELEASE FT33F001 

and then 

DDEF FT33F001 , VS, BBB 


and then type “GO” to continue execution. 
Subsequently, the user can request function 202, which 
will read in the desired data from dataset BBB. 

Series 200 - Data Input and Revision 

This series of nine functions is used for inputting basic 
problem-defining data from datasets or for inputting 
changes to that data. The changes are typed in from the 
user’s terminal. 

201 -Forming matrices for test case I. - Function 201 
is of use mainly when executing test case I, which is 
presented in appendix C. Function 201 simply forms all 
matrices and defines all dimensions, all of which are 
included in NAMELIST MATDAT. The specific 
matrices that this function forms for the test case are 
shown in table V. The dimensions for this problem are 
N = 3, NC = 2, NM= 1, NO-2, and ND = 2. 

202 -Primary function for problem-defining data 
input. - Function 202 is the one usually used for reading 
in (from a VS dataset) the data that define the user’s 
problem. The data must reside in the dataset in the 
following NAMELIST form: 

&MATDAT (data for variables in NAMELIST 

MATDAT -see table II) 

&END 

&REFS (data for variables in NAMELIST REFS - 

see table II) 

&END 

Note that data for the two NAMELIST’s must be in the 
order shown. It is not necessary to have any of the REFS 
parameters entered if the user wishes to have AESOP use 
the default values shown in table III. However, the entry 
in the aforementioned dataset must then appear as 

&REFS &END 

203, 204 , and 210— Functions used for revising data 
that define the user’s problem. - Functions 203, 204, and 
210 are useful when conducting design iterations, for 
filter design, regulator design, or plant sensitivity studies. 
Each function allows various parameter changes to be 
made via NAMELIST’s from the user’s terminal. (For 
NAMELIST definitions, see table II.) 

Function 203 is used to revise elements in the control 
design performance index weighting matrices QC, NN, 
and PCINV. NAMELIST CONPAR is used for this 
purpose. Function 204 is used to revise elements in noise 
power spectral density matrices QQ and RRINV through 
NAMELIST ESTPAR. These noise parameters, of 
course, strongly influence the characteristics of the 
associated Kalman filter. Function 210 is useful when 
changes are to be made in any of the variables in 
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TABLE IV. - SUMMARY OF AESOP FUNCTIONS 


Function 

Description 


Function 

Description 


Function 

Description 

Program control 

Compute 

Plot 

Store 

Transient function zeros 

101 

Change reference values 



Open-loop system analysis 




Open-loop system: 

102 

Program pause 








701 

u to z 




401 



Eigenvalues 



702 

u to y 

Inputting or revising matrices 


402 



Eigenvectors and mode shapes 



703 

w to z 




403 



Controllability, observability, 



704 

w to y 


Matrix input: 





and residues 



705 

Optimal controller 

201 

Matrices for test case 


404 



Normalization 





202 

Basic matrix data input 


405 



Unnormal ization 



Function 

Description 


Revising matrices: 




Frequency responses 


LQR 

Kalman filter 













cU J 

xjl, Nil, ana rcxNV por liikj 










204 

QQ and RRINV (for filter) 





Open loop: 



LQR, Kalman filter, and covariances 

210 

A, B, C, DOUT, H, CSP, QC, 


501 

502 

503 

u to z 






NN, PCINV, QQ, and RRINV 


504 

505 

506 

u to y 


801 

809 

Solve Riccati equation 




507 

508 

509 

w to z 


808 

816 

Store solution matrix SS/PP 


Reading stored matrices: 


510 

511 

512 

w to y 


802 

810 

Store gain KC/KE 

205 

KC (from LQR) 







805 

813 

Check positive-definiteness of SS/PP 

206 

KE (from filter) 





With state feedback: 


806 

814 

Check symmetry of SS/PP 

207 

SS (from LQR) 


513 

514 


w to y 


807 

815 

Compute residual error 

208 

PP (from filter) 


515 

516 


w to u 





209 

KFF 





With Kalman filter feedback: 



Function 

Description 












Matrix formation 


517 

Cl Q 

518 


w to z 



Eigenvalues and eigenvectors 




oiy 

DcU 


w to y 





301 

A - B*KC 


521 

522 


w to u 



803 

E-values of A - B*KC 

302 

A - B-KC - KE*H 


523 

524 

525 

Optimal controller 



811 

E-values of A - B*KC - KE*H 

303 

"Total" system matrices 








812 

E-values of ATOT 






Transient response 



804 

E-vectors of A - B*KC 




601 



Open-loop step 



Covariances, etc. 




602 



Open-loop initial condition 








603 



Closed-loop LQR initial condition 



817 

Covariances of LQR system 




604 



Nonzero-set-point LQR step 



818 

Covariance error check 









819 

Feedforward gain calculation (KFF) 









820 

Store KFF 






i 

User-supplied subroutines 








901 

UZR901 








902 

UZR902 








903 

UZR903 





i 



904 

UZR904 


a 



TABLE V. - INPUT MATRICES FOR THIRD-ORDER TEST CASE 


A = 

Input matrix 

CSP = 

Input matrix 


1 

2 

3 


1 

2 

3 

1 

-0. 1000D-00 

1.000 

0.0000 

1 

1.000 

0.0000 

0.0000 

2 

0.0000 

0.0000 

1.000 

2 

0.0000 

0.0000 

2.000 

3 

0.0000 

- 1.000 

-0.2000D-07 





B = 

1 

2 


QC = 

1 

2 

3 

1 

0.0000 

0.0000 


1 

500.0 

0.0000 

0.0000 

2 

0.0000 

1,000 


2 

0.0000 

9.000 

0.0000 

3 

1.000 

0.0000 


3 

0.0000 

0.0000 

0. 4000D-07 

D = 

1 

2 


HH = 

1 

2 


1 

0.0000 

0.0000 


1 

40.00 

0.0000 


2 

0.0000 

1.000 


2 

0.0000 

81.00 


3 

1.000 

0.0000 


3 

32.00 

0.0000 


C = 

1 

2 

3 

PCINV = 

1 

2 


1 

1.000 

0.0000 

0.0000 

1 

0. 1563D-01 

0.0000 


2 

0.0000 

0.0000 

1.000 

2 

0.0000 

0. 1235D-02 


H =: 

1 

2 

3 

QQ - 

1 

2 

3 

1 

1.000 

0.0000 

0.0000 

1 

0.0000 

0.0000 

0.0000 





2 

0.0000 

2.000 

0.0000 





3 

0.0000 

0.0000 

20.00 

DOUT = 

1 

2 


RRINV = 

1 



1 

0.0000 

0.0000 


1 

1.000 



2 

o.oooo 

1.000 







NAMELIST MATDAT, in particular, problem dimen- 
sions and matrices A, B, C, D, DOUT, H, and CSP. 
However, the noise matrices and the weighting matrices 
can be modified here also if so desired. 

205, 206, 207, 208, and 209 — Functions for reading 
gain and Riccati solution matrices. — AESOP contains 
three functions (801, 809, and 819) that solve Riccati 
equations and associated gain matrices and five functions 
(802, 808, 810, 816, and 820) that store results in datasets. 
Functions 205 to 209 are used for reading in datasets that 
contain gain or Riccati solution matrices previously 
computed by function 801, 809, or 819. Function 205 
reads control gain matrix KC from dataset CG$1 ($1 is 
the parameter in PROCDEF AESRUN). Function 206 
reads Kalman filter gain matrix KE from dataset EG$1. 
Function 207 reads control Riccati solution matrix SS 
from dataset SS$1. Function 208 reads Kalman filter 
Riccati solution matrix PP from dataset PP$1. Function 
209 reads feedforward gain matrix KFF from dataset 
FFG$1 . Note that by using PAUSE function 102, the user 
can appropriately re-datadef any of the preceding 
datasets so as to read in the data from the particular 
dataset desired. Refer to table I for definition of the unit 
numbers associated with these five datasets. 

Series 300 -Matrix Formation 

The three functions in series 300 all take gain and open- 
loop plant matrices and form various matrices related to 


closed-loop control system configurations. Separate 
functions are assigned to forming these matrices because 
a number of AESOP functions require these same 
matrices as input. Table VI summarizes the three 
functions and indicates for which other AESOP function 
these three functions are prerequisites. Examples of using 
these functions are given in appendix C; for 301 , pages 69 
and 79; for 302 and 303, page 70. 

Series 400 - Open-Loop System Analysis 

This series of five functions is used in conducting pre- 
design analyses for the open-loop plant or for plant data 
normalization. 

401 — Open-loop system eigenvalues.- Function 401 
simply computes the eigenvalues of the A matrix. The 
eigenvalues are displayed in both Cartesian (a±j(3) form 
and polar (frequency and f) where the frequency is in 
hertz and f is the damping ratio, f 4-cos(arctan |j3|/a). 
A negative f indicates a right-half-plane eigenvalue (or 
pair). An example appears on page 68 of appendix C. 

402 — Open-loop eigenvectors and mode 
shapes. - Function 402 computes and prints out the 
modified eigenvector matrix for the A matrix and also 
these vectors in mode-shape form. Refer to the section 
Theoretical Background for a description of the modified 
eigenvector and mode-shape format. Note that one must 
call function 401 before calling function 402. Use of 
function 402 is also shown on page 68 of appendix C. 
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TABLE VI. - AESOP FUNCTIONS USED FOR MATRIX FORMATION 


Function 

Matrix formed 

Dimension 

Equation 

Function where used 

301 

AMBKC 

NxN 

AMBKC « A - B-KC 

513,515,603, 

604,803,819 

302 

ABKCEH 

NxN 

ABKCEH = A - B-KC - KE-H 

523,705,811 

303 

AT0T, DTOT, CTOT, 
KCTOT, HTOT 

(a) 

(a) 

517,519,521,812 


a ATOT = j^ rR | ~KE r H J ’ 2Nx2N 

DTOT = [-§-]» 2NxND 

CTOT = [C i -DOUT'KC], NOx2N 
HTOT = [H i 0], NMx2N 
KCTOT = [0 i -KC], NCx2N 

403 - Controllability, observability, and 
residues . - Prerequisites for function 403 are functions 
401 and 402. Given the open-loop system described by A, 
B, C, and H, function 403 computes the controls effec- 
tiveness matrix (T -1 B) and two system observability 
matrices (HT and H C). Also, this function computes 
the residues for systems (A, B, C) and (A, B, H). All 
outputs are generated in mode-shape format. See page 68 
of appendix C for an example of the use of function 403. 

404— Normalization of system matrices. - Function 
404 allows the user to normalize all system matrices by 
using scaling factors that have been prestored in a 
dataset. The program prompts the user to datadef to unit 
34 the dataset containing these (normalizing) factors. The 
factors are defined as 


System 

Vector of 

Size 

variable 

normalization 


(vector) 

factor 


X 

sex 

50 

u 

see 

5 

z 

sez 

5 

y 

SCY 

50 

y.p 

SCYSP 

5 


The program initializes all normalizing factors to unity 
before reading in the dataset. The (VS) dataset must 
contain the data in NAMELIST form, where the 
NAMELIST name is NRMS. An example of a dataset’s 
contents might be 

&NRMS SCX(1) = 2., 3.,SCU(1) = 42., SCU(3) = .22 
&END 

Here only the normalizing factors for the first and second 
state variables and the first and third control variables are 
to be nonunity. Matrices affected by normalization are 


A, B, C, H, QQ, RRINV, D, DOUT, and CSP. 

Examples of using function 404 and companion function 
405 appear on pages 77 and 79 of appendix C. 

405— Unnormalization of matrices . - Function 405 is a 
companion to function 404 and allows resultant matrices 
KC, KE, KFF, and PP to be transformed back to 
dimensional form if desired. If prior normalization of 
input matrices is performed (via 404) in the same run of 
AESOP, no new scaling factors need be read in. 
However, function 405 automatically prompts the user 
(as is done in function 404) for the dataset containing the 
normalizing factors if not previously read in. 

Series 500 -Frequency Responses and Bode Plots; and 
Series 700— Transfer Functions 

Functions in these two series (500 and 700) are closely 
related in that all allow the user to perform frequency 
domain input/output calculations for both open- and 
closed-loop system configurations. A general description 
of these calculations is provided in the section Theoretical 
Background under the heading Transfer Functions and 
Frequency Response Calculations. All functions in series 
500 and 700 are summarized in table VII. Four system 
configurations are considered here: 


Configuration 

Figure 

Open loop 

8 

State-variable feedback 

9 

Kalman filter feedback 

10 

Optimal controller only 

11 


The figures describing the configurations appear in the 
section Theoretical Background. As noted in table VII, 
the user specifies the component of the input/output pair 
of interest by selecting values for indices JINC, JIND, 
IMEAS, or IOUT, which are all members of NAMELIST 
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TABLE VII. - AESOP FUNCTIONS FOR COMPUTING TRANSFER FUNCTIONS AND FREQUENCY RESPONSES 


System 
configuration 


Input 

variable 


Output 

variable 


AESOP program function numbers 





— 



Transfer functions 

Frequency responses 


Name 

Index 3 

Name 

Index 3 

— 

— 

. — 

__ — 








Numerator 

Poles 

Magnitude, phase. 

Response 

Store magnitude 







zeros and gain 


and coeff icientsb 

plots 

and phase 

Open loop 

u 

JINC 

z 

IMEAS 

701 

401 

501 

502 

503 

i 

u 

JINC 

y 

I0UT 

702 

I 

504 

505 

506 

1 

w 

JIND 

z 

IMEAS 

703 

1 

507 

508 

509 

T 





y 

I0UT 

704 

♦ 

510 

511 

512 

State-variable feedback 





y 

I0UT 

— 

803 

513 

c 514 

— 

State-variable feedback 





u 

JINC 



803 

515 

516 

— 

Kalman filter feedback 





z 

IMEAS 



812 

517 

c 518 

— 

Kalman filter feedback 





y 

I0UT 

— 

812 

519 

c 520 

— 

Kalman filter feedback 





u 

JINC 



812 

521 

522 

— 

Optimal controller only 

z 

IMEAS 

u 

JINC 

705 

811 

523 

524 

525 


a Indices are included in NAMELIST REFS (table III). Index values determine vector input/output pair for which the 
transfer function or frequency response will be computed. 

Computes (1) coefficients of the transfer function numerator and denominator polynomials and (2) frequency response 
magnitudes and phase angles for initial frequency FI(Hz), frequency spacing DELF(Hz), number of points IF(<JL000), and 
a time scale factor TSFTR, 

c If NCURV = 2, this closed-loop response will be crossplotted with its corresponding open-loop response. Note that 
the open-loop response used is the last one calculated. 


REFS. It can be seen that transfer function numerator 
and denominator polynomial coefficients (eq. (42)), 
frequency response calculations, storing of frequency 
responses, and Bode plotting are available for all four 
system configurations. Specific function numbers allow 
the user to choose between control u or disturbance w 
inputs and between noisy z or noise-free y outputs. One 
exception is for the optimal-controller-only configura- 
tion, where the only input is measurement vector z and 
the output is the control vector u. Examples of the use of 
500 series functions appear on pages 72 to 74 of appendix 
C, including plotted output. Transfer function zeroes and 
gain information can be obtained only for the open-loop 
plant and for optimal controller configurations. How- 
ever, transfer function poles can be computed for all 
configurations. Note that these are all eigenvalue 
computations and are included in series 400 or series 800. 
Examples of the use of functions 701 to 705 are shown on 
pages 76 and 77 of appendix C. 

Frequency response data (frequency, amplitude, and 
phase angle), which are stored on datasets by functions 
503, 506, 509, 512, and 525, are all written by using a 
Fortran WRITE statement of the form 

WRITE (i, b) (FREQ(I),AMP(I), PHASE(I),I = I,IF) 

Here, i is an appropriate unit number, IF is the number 
of response data points, and FORMAT statement b is 
FORMAT (10G12.5). An appropriately formatted 
READ statement would be used if one wished to use any 
of these computed frequency responses in a separate 
program. 


Series 600 - Transient Responses 

Transient responses - for either step or initial 
condition inputs -are computed by the four AESOP 
functions in this series. Pertinent information required to 
use these functions is shown in table VIII. Functions 601 
and 602 are for the open-loop system configuration and 
603 and 604, for linear quadratic regulators. Note that to 
use a function, the user should specify 

(1) Input amplitude (or amplitudes), AMP 

vectors 

(2) Desired time step, DT 

(3) The number of time points to be computed, 
ITRMX < 1000 

(4) Desired input/output response pairs for which 

responses are to be calculated, MS or MIC 

All of these parameters are in NAMELIST REFS and 
are set to default values in the main AESOP program (see 
table III for default values). The transient response plots 
that are generated by these functions will be generated on 
the particular graphics device that the user has previously 
defined. Examples of these plots and the use of functions 
601 to 604 appear in test case I in appendix C. 

Series 800 - LQR and Filter Design 

Functions in series 800 are for (1) solving the Riccati 
equations associated with the LQR or Kalman filter 
design problems or (2) related gain, covariance, or 
eigenvalue/eigenvector calculations. 

Functions for solving Riccati equations. -Table IX 
lists the numbers of the functions in the 800 series that 
deal directly with LQR or Kalman filter design. Three 
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TABLE VIII. - AESOP FUNCTIONS FOR COMPUTING TRANSIENT RESPONSES 


Function 

Type of transient 

Input 

Vector containing 

Output 

Input/output 



variable 

input amplitudes 

variables 

select matrix 3 

601 

Step response of 

u 

AMPSR 

y 

MSR0LY 


open-loop system 



X 

MSR0LX 

602 

Initial-condition 



y 

MICOLY 


response of open- 

x(0) 

AMPICX 




loop system 



X 

MIC0LX 

603 

Initial-condition 



y 

MICCLY 


response of linear 

x(0) 

AMPICX 

X 

MICCLX 


quadratic regulator 



u 

MICCLU 

604 

Step response of 



y$p 

MSPYSP 


nonzero-set-point 

•Yspd 

AMPSP 

y 

MSPY 


regulator 


u 

MSPU 


a Use of input/output select matrices: 


M XXXX (input index, output index) = 1, if response is to be calculated and plotted 

0, otherwise 

That is, for function 601, if response of y(3) to u(2) is desired, set MSR0LY(2,3) = 1 
or, for function 602, if response of y(l) to an initial condition on x(2) is desired, 
set MICOLY (2,1) =1. All amplitude vectors and input/output select matrices are de- 
faulted to all ones. They can be changed via NAMELIST REFS. All responses are for 
only one input or initial-condition component applied. Default for time step DT is 
0.01 and that for time points ITRMX is 100. They can be changed via NAMELIST REFS. 


TABLE IX. - RICCATI EQUATION SOLUTION 


Task 

LQR design 
problem 

Kalman filter 
design problem 


Function 

Solve matrix Riccati equation: 



Number of function 

801 

809 

Name of solution matrix 

SS 

PP 

Name of gain matrix 

KC 

KE 

Store Riccati results in dataset: 



Function for storing solution matrix 

808 

816 

Function for storing gain matrix 

802 

810 

Riccati equation accuracy checks 
performed on solution matrix: 



Positive-definiteness 

805 

813 

Symmetry 

806 

814 

Residual error matrix 

807 

i 

815 


types of functions are provided: (1) those for solving 
Riccati equations, (2) those for storing Riccati equation 
results in datasets for future use, and (3) those for 
checking the accuracy of Riccati solutions. Companion 
functions are provided in series 200 for subsequently 
reading in the matrices computed and stored by functions 
in the 800 series. Examples of the use of functions 
involved with the LQR problem (801 to 808) are given in 
appendix C, pages 68 to 70. Examples of Kalman filter 
design calculations (functions 809 to 816) appear on 
pages 70 to 72 of appendix C. 

To read Riccati or gain matrices from datasets into 
another program (e.g., into a program that implements 
an LQR control law), the user must know the correct 
Fortran READ statements to use. The appropriate 
(unformatted) READ statements are as follows: 


(1) For reading the control Riccati solution (stored in 
dataset SS$1): 

READ (i)((SS(I, J),I = 1 ,N),J = 1 ,N) 

(2) For reading the Kalman filter Riccati solution 
(stored in dataset PP$1): 

READ (i)((PP(I,J),I= 1,N),J= 1,N) 

(3) For reading the LQR gains (stored in dataset 
CG$1): 

READ (i)((KC(I,J),I = 1 ,NC),J = 1 ,N) 

(4) For reading the Kalman filter gains (stored in 
dataset EG$1): 
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READ (i)((KE(I,J),I = 1,N),J = 1,NM) 

(5) For reading the feedforward gains (stored in 
dataset FFG$1): 

READ (i)((KFF(I,J),I = 1,NC),J = 1,NC) 

It is assumed that the user datadefs unit “/” to the 
appropriate dataset. 

Functions for eigenvalue /eigenvector compu- 
tation. -Table X lists the functions that compute eigen- 
values and eigenvectors associated with either regulator 
or filter designs. The eigenvalues and eigenvectors of the 
matrix A-BKC are those of the linear regulator. The 
eigenvalues of matrix A-BKC-KEH are those of the 
optimal controller (depicted in fig. 11). The eigenvalues 
of matrix ATOT are those of A-BKC plus those of 
A - KEH, the latter being the Kalman filter eigenvalues. 


TABLE X. - EIGENVALUE AND EIGENVECTOR 
COMPUTATIONS IN SERIES 800 


Matrix 

Eigenvalue 

Eigenvector 


Function 

A - B*KC 

803 

804 

A - B*KC - KE*H 

811 

— 

ATOT a 

812 

— 

a Matrix defined 

in table VI. 



817 -Covariance matrices for system controlled by 
linear stochastic regulator. - Function 817 implements 
the solution to the state covariance matrix (Lyapunov) 
equation given by equation (32) plus associated matrix 
equations (33), (34), and (35). The latter equations 
compute control, measurement, and output covariance 
matrices. Normally, function 817 is called after com- 
puting LQR gain matrix KC and Kalman filter gain 
matrix KE. However, one or both of these two matrices 
may be zero when function 817 is called. The resultant 
state covariance matrix XX will then correspond to the 
meaningful system configurations as shown in the follow- 
ing table. An example of the use of function 817 appears 
on page 72 of appendix C. 


LQR gain 
matrix, 

KC 

Kalman gain 
matrix, 

KE 

System configuration for which 
XX is covariance matrix 

Figure 

*> 

*) 

Linear stochastic regulator 

7 

= 0 

= 0 

Open -loop system with noise input 

3 

/o 

-0 

State feedback system with 
noise input 

9 

= 0 

/O 

Open-loop system with Kalman 
filter 

6 


818— Covariance matrix error check. - Function 818 
may be called after function number 817 to compute an 
estimate of the error incurred in computing the state 
covariance matrix. An example of its use appears on page 
72 of appendix C. 

819 — Feedforward gain matrix for nonzero-set-point 
regulator. - Function 819 computes feedforward gain 
matrix KFF, which is part of the nonzero-set-point 
regulator. Matrix KFF is given by 

KFF = [ - CSP(A - BKC) “ 1 B] -1 

The matrix KFF is needed when computing the closed- 
loop system step responses with function 604. Functions 
819 and 820 are demonstrated in test case I, appendix C, 
on page 72. 

820— Store feedforward gain matrix.- Function 820 
stores KFF in a dataset by using the following un- 
formatted WRITE statement: 

WRITE(iii)((KFF (I , J) , I = 1 ,NC), J = 1 ,NC) 

The user can read these data into another program by 
using a matching Fortran READ statement. 

Series 900 - User-Supplied Subroutines 

A series of four function numbers have been set aside 
for use as user-defined subroutines. These subroutines 
are called UZR901, UZR902, UZR903, and UZR904. 
Thus users may write their own special-purpose sub- 
routines by using any of the above as subroutine names. 
Linkage between the subroutine and the main AESOP 
program would be achieved by using any of the 
COMMON’S used in the AESOP program. Users need 
not recompile the AESOP subroutine containing CALLS 
to the four user-supplied subroutines, but need only 
compile their subroutines into a library that has higher 
priority in the JOBLIB chain than the ‘"standard” 
AESOP program library. 


Concluding Remarks 

The program described in this report was not meant to 
be a static entity but rather was envisioned as a basic 
structure to which can be added new and useful system 
design functions as the need arises. Some functions 
presently contemplated for future inclusion are discrete 
LQR and Kalman filter gain calculations, time responses 
for Kalman filters to noise signal inputs, linear model- 
order reduction procedures, LQR or Kalman filter eigen- 
value sensitivity calculations, transient responses of 
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discrete LQR or Kalman filters, and multivariable fre- 
quency domain control design algorithms. 

Modifications are also envisioned that could improve 
program-user interfacing. One would be to replace the 
present function-number-input system with a set of 
mnemonics (i.e., the user would enter LQR instead of 801 
when requesting an LQR problem solution, etc.). 
Another addition would be to allow the user to use a light 
pen to select AESOP functions from a menu displayed on 


a CRT terminal screen. The present program structure 
has sufficient flexibility to accommodate these types of 
additions without compromising present interactive 
capability. 


Lewis Research Center 

National Aeronautics and Space Administration 
Cleveland, Ohio, June 14, 1983 
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Appendix A 
Symbols 


Dimensions are given for all vectors, matrices, and 
higher dimensional arrays. Where dimensions are not 
given, the variable is a scalar quantity. 


Variable 

Dimension 

Description 

A 

NxN 

system matrix 

A* 

nxn 

matrix whose eigenvalues 
are transfer function 
zeroes 

A 

nxn 

general matrix 

ABKCEH 

NxN 

matrix A-B KC-KE H 

AMBKC 

N X N 

matrix A-BKC 

AMPICX 

N 

vector of initial-condition 
amplitudes 

AMPSP 

NC 

vector of input step 
amplitudes for closed-loop 
system 

AMPSR 

NC 

vector of input step 
amplitudes for open-loop 
system 

ATOT 

NTOT x NTOT system matrix for combined 
regulator-Kalman filter 
system 

' a k 


coefficient of transfer 
function numerator 
polynomial 

B 

NxNC 

control input matrix 

B 

N x NC 

controllability matrix 

B 

nxnc 

general matrix 

b k 


coefficient of transfer 
function denominator 
polynomial 

C 

NO X N 

output matrix 

C 

noxn 

general matrix 

CSP 

NCxN 

set -point output matrix 

CSP 

NCxN 

normalized CSP matrix 

CTOT 

NO x NTOT 

output matrix for combined 
regulator-Kalman filter 
system 

D 

noxnc 

general matrix 

D 

NxND 

disturbance input matrix 

DELF 


spacing between frequency 
points 

DOUT 

NOxNC 

feedforward matrix 

DT 


time step 


DTOT 

NTOT x ND 

disturbance input matrix 
for combined regulator- 
Kalman filter system 

e 

N 

Kalman filter estimation 
error vector 

E 

NxN 

error matrix (X - X) 

% 

NMxNC 

selection matrix 

FI 


initial frequency 

G(s) 

NMxNC 

transfer function matrix 

g(s) 


transfer function 

H 

NMxN 

measurement matrix 

H 

NMxN 

observability matrix 

HTOT 

NM X NTOT 

measurement matrix for 
combined regulator- 
Kalman filter system 

I 

NxN 

identity matrix 

IF 


integer, number of desired 
frequency response points 

IMEAS 


integer, index of meas- 
urement 

IOUT 


integer, index of output 

ISPACE 


integer, controls frequency 
of frequency response 
printouts 

ITRMX 


integer, number of desired 
time response points 

J 


performance index 

JINC 


integer, index of control 

JIND 


integer, index of distur- 
bance 



transfer function gain 

KC 

NCxN 

control gain matrix 

KCTOT 

NC X NTOT 

control gain matrix for 
combined regulator- 
Kalman filter system 

KE 

NxNM 

Kalman filter gain matrix 

KFF 

NCxNC 

feedforward gain matrix for 
nonzero-set-point 
regulator 

k* 


feedback gain term 

LINLOG 


integer, indicates whether 
frequency response plots 
are to be linear, log, or 
both 

e 


matrix dimension 
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MICCLU 

N x NC 

integer, input/output, 
closed-loop, initial- 
condition response 
selection matrix for 
controls; state feedback 

MICCLX 

NX N 

integer, input /output, 
closed-loop, initial- 
condition response 
selection matrix for 
states; state feedback 

MICCLY 

NxNO 

integer, input /output, 
closed-loop, initial- 
condition response 
selection matrix for 
outputs; state feedback 

MICOLX 

NX N 

integer, input /out put, 
open-loop, initial- 
condition response 
selection matrix for 
states 

MICOLY 

NxNO 

integer, input/output, 
open-loop, initial- 
condition response 
selection matrix for 
outputs 

MSPU 

NCxNC 

integer, input /out put, 
closed-loop, step response 
selection matrix for 
controls; state feedback 

MSPY 

NCxNO 

integer, input/output, 
closed-loop, step response 
selection matrix for out- 
puts; state feedback 

MSPYSP 

NCxNC 

integer, input /output, 
closed-loop, step response 
selection matrix for set- 
point outputs; state 
feedback 

MSROLX 

NCxN 

integer, input/output, 
open-loop, step response 
selection matrix for states 

MSROLY 

NCxNO 

integer, input/output, 
open-loop, step response 
selection matrix for 
outputs 

M 

Nxf 

matrix, general factor of 
QC matrix 

N 


integer, number of states 

NC 


integer, number of controls 

NCURV 


integer, controls cross 


plotting of frequency 
responses 


ND 


integer, number of 
disturbances 

NM 


integer, number of 
measurements 

NMAX 


integer, maximum number 
of states 

NN 

NxNC 

state-control weighting 
matrix 

NO 


integer, number of outputs 

NTOT 


integer, 2 x N 

n 


number of states for general 
dynamic system 

nc 


number of inputs for 
general dynamic system 

no 


number of outputs for 
general dynamic system 

PCINV 

NCxNC 

inverse of control weighting 
matrix 

PP 

NxN 

Kalman filter error 
covariance matrix 

Pk 


transfer function pole 

Q 

txt 

symmetric, positive-definite 
matrix 

QC 

NxN 

state weighting matrix 

QQ 

NXN 

power spectral density 
matrix of plant 
disturbance 


NMxNC 

residue matrix of transfer 
function matrix G(s) 

R 

NxN 

residual matrix 

RRINV 

NMxNM 

inverse of power spectral 
density matrix of 
measurement noise 

r J 


jth residue of transfer 
function g(s) 

scu 

NC 

vector of normalization 
factors for controls 

sex 

N 

vector of normalization 
factors for states 

SCY 

NO 

vector of normalization 
factors for outputs 

SCYSP 

NC 

vector of normalization 
factors for set-point 
outputs 

SCZ 

NM 

vector of normalization 
factors for measurements 

ss 

NxN 

Riccati solution matrix for 
linear quadratic regulator 

s 


Laplace variable 
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T 

NxN 

modified eigenvector 
matrix 

TSFTR 


time scale factor 

t 


time 

+ 1 

N 

modified eigenvectors 

uu 

NCxNC 

control covariance matrix 

u 

NC 

control vector 

v 

NM 

measurement noise vector 

W 

NxN 

general matrix 

w 

ND 

disturbance vector 

X 

NxN 

general matrix 

X 

NxN 

computed value of matrix X 

XX 

NxN 

state covariance matrix 

X 

N 

state vector 

X 

N 

modal state vector 

X 

N 

estimated state vector 

YY 

NOxNO 

output covariance matrix 

y 

NO 

output vector 

y SP 

NC 

set-point vector 

y sP d 

NC 

desired set -point vector 

ZERMAX 


maximum expected value of 
transfer function zeroes 

zz 

NMxNM 

measurement covariance 
matrix 

z 

NM 

measurement vector 

Zl 

NM 

noise-free component of 
vector z 


z NM 

Zk 

a 

0 

r NxNC 

7 N 

<5 N 

<5(r) 

f 

7 7 N 

A NxN 

\»\*+ 1 
r 

* NxN 

Superscripts: 

T 

-1 

Operators: 

E[ } 


estimated value of z 
transfer function 
numerator zero 
real part of eigenvalue 
imaginary part of 
eigenvalue 

forced-response matrix for 
discrete system 
real part of eigenvector 
imaginary part of 
eigenvector 
unit impulse function 
damping ratio 
eigenvector 

block diagonal form of A 
matrix 

complex eigenvalue pair 
time 

state-transition matrix 

matrix transpose 
matrix inverse 

expected value of 
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Appendix B 

Subroutine Descriptions 


This appendix describes the AESOP main program and 
all associated subroutines. A short description of each 
subroutine is given, including what subroutines may call 
it, and what subroutines it calls. Where a subroutine has 
a parameter list, all variables in that list are defined. 

A tape of the AESOP program documented in this 
report is available from COSMIC. That version of 
AESOP has been sized to accommodate systems having 
the following dimensions: 

N<50 
NM<5 
NC <5 
ND < 15 
NO <50 


If these dimensions are greater than or equal to the 
corresponding dimensions of the user’s problem, the user 
need make no changes to the main AESOP program 
(assuming that the program will fit on the user’s 
computer). However, the user may wish to resize the 
program, either to accommodate problems with larger 
dimensions or to reduce storage requirements in order to 
allow the program to fit on a smaller computer. As 
currently dimensioned, AESOP requires approximately 
300 000 bytes for all Fortran source codes and about 
1 700 000 bytes for variable storage (the variables in 
COMMON’S appearing in the main program). To change 
program dimensions, array dimensions appearing in 11 
labeled COMMON’S must be changed. These 
COMMON’S appear only in the main program and the 
nine main subroutines (AES 100 to AES900). 


AESOP 

AESOP IS THE MAIN ROUTINE. AESOP CALLS SUBROUTINES AES100, 

AES200, AES300, AES400, AES500, AES600, AES700, AES800, 

AES900, AND THE PLOTTING SUBROUTINES. 

AESOP IS NOT CALLED BY ANY SUBROUTINES. 

****************************************************************** 


SUBROUTINE AES100 (IFN, IFUNC, I AND, MZ, IPRT, WHEN) 
****************************************************************** 

SUBROUTINE AES100 CONTAINS THOSE FUNCTIONS WHICH PERFORM PROGRAM 
CONTROL. AES100 CALLS SUBROUTINE PREREQ. 

AES100 IS CALLED BY MAIN PROGRAM AESOP. 

INPUTS: 

IFN VECTOR OF FUNCTION NUMBERS TO BE DONE (1000) 

IFUNC FUNCTION NUMBER 

MZ WHICH FUNCTION IS TO BE DONE 

IPRT PRINT OPTION: 1, STANDARD PRINT; 

2, EXTENDED PRINT 

WHEN LOGICAL MATRIX OF PREREQUISITES (450,50) 

OUTPUTS: 

I AND DECISION VARIABLE: 

0, PREREQUISITES HAVE BEEN DONE; 

1, PREREQUISITES HAVE NOT BEEN DONE 



SUBROUTINE AES200 (IFN, IFUNC, I AND, MZ, IPRT, WHEN) 

******* ****************************** ***************************** 


SUBROUTINE A ES200 CONTAINS THOSE FUNCTIONS WHICH INPUT DATA. 

AES200 CALLS SUBROUTINES MATCHG, MATIN, MATPRT, MATRD, AND PREREQ. 
AES200 IS CALLED BY MAIN PROGRAM AESOP. 


INPUTS: 


IFN VECTOR OF FUNCTION NUMBERS TO BE DONE (1000) 
IFUNC FUNCTION NUMBER 

MZ WHICH FUNCTION IS TO BE DONE 

IPRT PRINT OPTION: 1, STANDARD PRINT; 

2, EXTENDED PRINT 

WHEN LOGICAL MATRIX OF PREREQUISITES (450,50) 

OUTPUTS: 

IAND DECISION VARIABLE: 

0, PREREQUISITES HAVE BEEN DONE; 

1, PREREQUISITES HAVE NOT BEEN DONE 


****************************************************************** 


SUBROUTINE AES300 (IFN, IFUNC, IAND, MZ, IPRT, WHEN) 
****************************************************************** 

SUBROUTINE AES300 CONTAINS THOSE FUNCTIONS WHICH FORM MATRICES. 
AES300 CALLS SUBROUTINES MATPRT AND PREREQ. 

AES300 IS CALLED BY MAIN PROGRAM AESOP. 

INPUTS: 

IFN VECTOR OF FUNCTION NUMBERS TO BE DONE (1000) 

IFUNC FUNCTION NUMBER 

MZ WHICH FUNCTION IS TO BE DONE 

IPRT PRINT OPTION: 1, STANDARD PRINT; 

2, EXTENDED PRINT 

WHEN LOGICAL MATRIX OF PREREQUISITES (450,50) 

OUTPUTS : 

IAND DECISION VARIABLE: 

0, PREREQUISITES HAVE BEEN DONE; 

1, PREREQUISITES HAVE NOT BEEN DONE 

****************************************************************** 


SUBROUTINE AES400 (IFN, IFUNC, IAND, MZ, IPRT, WHEN) 
****************************************************************** 

SUBROUTINE AES400 CONTAINS THOSE FUNCTIONS WHICH CALCULATE 
CONTROLLABILITY, OBSERVABILITY, EIGENVALUES, EIGENVECTORS, AND 
RESIDUES, AND WHICH PERFORM NORMALIZATION AND UN-NORMALIZATION. 
AES400 CALLS SUBROUTINES CTBL, EGVCTR, EIGEN, MATPRT, MODSHP, 

NRML, OBSBL, PREREQ, RESI, AND UNRML. 

AES400 IS CALLED BY MAIN PROGRAM AESOP. 

INPUTS: 


IFN 

IFUNC 

MZ 


VECTOR OF FUNCTION NUMBERS TO BE DONE (1000) 

FUNCTION NUMBER 

WHICH FUNCTION IS TO BE DONE 



IPRT PRINT OPTION: 1, STANDARD PRINT; 

2, EXTENDED PRINT 

WHEN LOGICAL MATRIX OF PREREQUISITES (450,50) 

OUTPUTS: 

IAND DECISION VARIABLE: 

0, PREREQUISITES HAVE BEEN DONE; 

1, PREREQUISITES HAVE NOT BEEN DONE 

****************************************** ************** ********* * 


SUBROUTINE AES500 (IFN, IFUNC, IAND, MZ, IPRT, WHEN) 
****************************************************************** 

SUBROUTINE AES500 CONTAINS THOSE FUNCTIONS WHICH CALCULATE 
FREQUENCY RESPONSE AND BODE PLOTS. 

AES500 CALLS SUBROUTINES BODE, FRSPNS, AND PREREQ. 

AES500 IS CALLED BY MAIN PROGRAM AESOP. 

INPUTS: 


IFN VECTOR OF FUNCTION NUMBERS TO BE DONE (1000) 

IFUNC FUNCTION NUMBER 

MZ WHICH FUNCTION IS TO BE DONE 

IPRT PRINT OPTION: 1, STANDARD PRINT; 

2, EXTENDED PRINT 

WHEN LOGICAL MATRIX OF PREREQUISITES (450,50) 

OUTPUTS: 

IAND DECISION VARIABLE: 

0, PREREQUISITES HAVE BEEN DONE; 

1 , PREREQUISITES HAVE NOT BEEN DONE 

****************************************************************** 


SUBROUTINE AES600 (IFN, IFUNC, IAND, MZ, IPRT, WHEN) 
****************************************************************** 

SUBROUTINE AES600 CONTAINS THOSE FUNCTIONS WHICH CALCULATE 
TIME RESPONSES AND THE ASSOCIATED PLOTS. 

AES600 CALLS SUBROUTINES DSCRT, ICRSP, MATPRT, PREREQ, AND STP. 
AES600 IS CALLED BY MAIN PROGRAM AESOP. 

INPUTS: 

IFN VECTOR OF FUNCTION NUMBERS TO BE DONE (1000) 

IFUNC FUNCTION NUMBER 

MZ WHICH FUNCTION IS TO BE DONE 

IPRT PRINT OPTION: 1, STANDARD PRINT; 

2, EXTENDED PRINT 

WHEN LOGICAL MATRIX OF PREREQUISITES (450,50) 

OUTPUTS: 

IAND DECISION VARIABLE: 

0, PREREQUISITES HAVE BEEN DONE; 

1 , PREREQUISITES HAVE NOT BEEN DONE 


★★★★★★****★★★* * ******** * ********* *** ****** A A ** * ******************* 



SUBROUTINE AES700 (IFN, IFUNC, I AND, MZ, IPRT, WHEN) 


SUBROUTINE AES700 CONTAINS THOSE FUNCTIONS WHICH CALCULATE 
ZEROES AND GAINS, 

AES700 CALLS SUBROUTINES GAIN, MATPRT, PREREQ, AND ZEROES. 
AES700 IS CALLED BY MAIN PROGRAM AESOP. 


INPUTS: 


IFN VECTOR OF FUNCTION NUMBERS TO BE DONE (1000) 
IFUNC FUNCTION NUMBER 
MZ WHICH FUNCTION IS TO BE DONE 

IPRT PRINT OPTION: 1, STANDARD PRINT; 

2, EXTENDED PRINT 

WHEN LOGICAL MATRIX OF PREREQUISITES (450,50) 

OUTPUTS: 


I AND DECISION VARIABLE: 

0, PREREQUISITES HAVE BEEN DONE; 

1, PREREQUISITES HAVE NOT BEEN DONE 


SUBROUTINE AES800 (IFN, IFUNC, I AND, MZ, IPRT, WHEN) 

SUBROUTINE AES800 CONTAINS THOSE FUNCTIONS WHICH CALCULATE 
EIGENVALUES AND EIGENVECTORS OF CONTROLS OR FILTERS, AND THE 
FEED FORWARD, RICCATI, AND COVARIANCE EQUATIONS. 

AES800 CALLS SUBROUTINES CONTRL, COVAR, EGVCTR, EIGEN, ESTMAT, 
LYPCK, MATPRT, MODSHP, MXINV, PREREQ, AND RICCHK 
AES800 IS CALLED BY MAIN PROGRAM AESOP. 

INPUTS: 

IFN VECTOR OF FUNCTION NUMBERS TO BE DONE (1000) 
IFUNC FUNCTION NUMBER 

MZ WHICH FUNCTION IS TO BE DONE 

IPRT PRINT OPTION: 1, STANDARD PRINT; 

2, EXTENDED PRINT 

WHEN LOGICAL MATRIX OF PREREQUISITES (450,50) 

OUTPUTS; 

I AND DECISION VARIABLE: 

0, PREREQUISITES HAVE BEEN DONE; 

1, PREREQUISITES HAVE NOT BEEN DONE 


SUBROUTINE AES900 (IFN, IFUNC, IAND, MZ, IPRT, WHEN) 
*★★**★★★★****★**★**★★*****★★★***★*★★****★★***★*★★★★★★****★★★*★★★★★ 

SUBROUTINE AES900 CONTAINS THOSE FUNCTIONS WHICH ARE SUPPLIED 
BY THE USER. AES900 DOES NOT CALL ANY SUBROUTINES. 

AES900 IS CALLED BY MAIN PROGRAM AESOP. 

INPUTS: 


IFN VECTOR OF FUNCTION NUMBERS TO BE DONE (1000) 
IFUNC FUNCTION NUMBER 



MZ WHICH FUNCTION IS TO BE DONE 

IPRT PRINT OPTION: 1, STANDARD PRINT; 

2, EXTENDED PRINT 

WHEN LOGICAL MATRIX OF PREREQUISITES (450,50) 

OUTPUTS: 


IAND DECISION VARIABLE: 

0, PREREQUISITES HAVE BEEN DONE; 

1 , PREREQUISITES HAVE NOT BEEN DONE 

****************************************************************** 


SUBROUTINE ARRAY ( IOPT, I , J,NROW,A) 
***************************************************************** * 

SUBROUTINE ARRAY CONVERTS AN ARRAY FROM VECTOR TO MATRIX OR THE 
REVERSE, ARRAY DOES NOT CALL ANY SUBROUTINES. ARRAY IS CALLED BY 
SUBROUTINES COVAR, EGVCTR, EIGEN, LYPCK, RICSS, AND ZEROES. 

INPUTS: 


IOPT OPTION INDICATING TYPE OF CONVERSION 

1 - FROM VECTOR TO MATRIX 

2 - FROM MATRIX TO VECTOR 

I NUMBER OF ROWS IN ACTUAL MATRIX 

J NUMBER OF COLUMNS IN ACTUAL MATRIX 

NROW NUMBER OF ROWS SPECIFIED FOR THE MATRIX A IN 
DIMENSION STATEMENT 

A IF MODE = 1, CONTAINS A VECTOR OF I*J LENGTH. 

IF MODE = 2, CONTAINS A MATRIX OF N BY J SIZE. 

OUTPUTS: 

A IF MODE = 1, CONTAINS A MATRIX OF N BY J SIZE. 

IF MODE = 2, CONTAINS A VECTOR OF I*J LENGTH. 

* ***************************************************************** 


BLOCK DATA 

**************************** ****************************** ******* 

THE BLOCK DATA SUBROUTINE CONTAINS ONLY INFORMATION FOR PLOT 
TITLES AND LABELS 

***************************************************************** 


SUBROUTINE BODE (FRQ, Al, A2, PHI1, PHI2, TTITL, TB, NPTS, KI, 

1 AMP, PHA, SETAP, KTYPE1, KTYPE2, IP, NAME, IONPLT) 

****************************************************************** 

SUBROUTINE BODE MAKES PLOTS OF FREQUENCY RESPONSES. AMP, PHA AND 
SETAP MUST BE SINGLE PRECISION BECAUSE OF THE PLOT SUBROUTINES. 
BODE CALLS PLOTTING SUBROUTINES ONLY. BODE IS CALLED BY 
SUBROUTINE AES500. 

INPUTS: 


VECTOR OF FREQUENCY 

(DIMENSION IS LESS THAN OR EQUAL TO 500) 
VECTOR OF AMPLITUDE FOR 1ST CURVE 
(DIMENSION IS LESS THAN OR EQUAL TO 500) 


FRQ 

Al 



A2 VECTOR OF AMPLITUDE FOR 2ND CURVE, IF DESIRED 
(DIMENSION IS LESS THAN OR EQUAL TO 500) 

PHI 1 VECTOR OF PHASE FOR 1ST CURVE 

(DIMENSION IS LESS THAN OR EQUAL TO 500) 

PH12 VECTOR OF PHASE FOR 2ND CURVE, IF DESIRED 
(DIMENSION IS LESS THAN OR EQUAL TO 500) 

TTITL TITLE OF PLOT 

(DIMENSION IS GREATER THAN OR EQUAL TO 15) 

TB TITLE OF PLOT 

(DIMENSION IS GREATER THAN OR EQUAL TO 15) 

NPTS NUMBER OF POINTS PER CURVE 
(LESS THAN OR EQUAL TO 500) 

KI 1, IF ONE CURVE PER PLOT 

2, IF TWO CURVES PER PLOT 
KTYPE1 ) THESE TWO VARIABLES DEFINE WHETHER 
) THE PLOT IS TO BE LINEAR, 

KTYPE2) LOG, OR SEMI-LOG 

IP PLOT ENTITY INDEX (USED BY PLOTSUBS ONLY) 
INCREASES BY ONE FOR EACH FRAME 
NAME NAME OF PLOT DATASET (9) (USED BY PLOTSUBS ONLY) 
(PARTITIONED DATASET THAT HOLDS PLOT ENTITIES) 
IONPLT 0, IF OFFLINE PLOTS 
1, IF ONLINE PLOTS 

OUTPUTS: 

AMP STORAGE VECTOR OF AMPLITUDES FOR 1 OR 2 CURVES 
(DIMENSION IS LESS THAN OR EQUAL TO 1000) 

PHA STORAGE VECTOR OF PHASES FOR 1 OR 2 CURVES 
(DIMENSION IS LESS THAN OR EQUAL TO 1000) 

IP PLOT ENTITY INDEX (USED BY PLOTSUBS ONLY) 
INCREASES BY ONE FOR EACH FRAME 

TEMPORARY STORAGE: 

SETAP VECTOR 

(DIMENSION IS LESS THAN OR EQUAL TO 500) 
****************************************************************** 


SUBROUTINE BOLLIN (A, AS, 8, C, X, Y, Zl, Z2, Z3, N, NMAX) 

★ ****★****★*★*★*****★★*★*****★*****★★***★***★*★**■★**★*★**★***★■*■•*■** 

SUBROUTINE BOLLIN CONVERTS X(DOT)=AX+BU, Y=C(TRANSPOSE)X TO 
TRANSFER FUNCTION Y/U=Z2/Z1 RATIO OF POLYNOMIALS. 

BOLLIN CALLS SUBROUTINES DAVISO AND DANSKY. BOLLIN IS CALLED BY 
SUBROUTINE FRSPNS. 

INPUTS: 

A SYSTEM MATRIX (N,N) 

B SYSTEM VECTOR (N) 

C SYSTEM VECTOR (N) 

N ACTUAL SIZE OF MATRIX A 

NMAX MAXIMUM SIZE OF N 

OUTPUTS: 

Zl DENOMINATOR COEFFICIENT VECTOR (N) 

12 NUMERATOR COEFFICIENT VECTOR (N) 

Z3 NUMERATOR COEFFICIENT VECTOR (N) 

TEMPORARY STORAGE: 

AS MATRIX (N,N) 

X VECTOR (N) 

Y VECTOR (N) 


****************************************************************** 



SUBROUTINE CONDI (VARO, SS, S, IN, JBL, IOR, NBL, IBL, IC, D 
1 IOPI, N, NMAX) 


SUBROUTINE CONDI CHANGES CONDITION OF A MATRIX BY PUTTING IT IN 
BLOCK DIAGONAL FORM (IF REDUCIBLE) AND THEN SCALING. 

CONDI CALLS SUBROUTINES REDU AND SCALEA. CONDI IS CALLED BY 
SUBROUTINES EIGEN AND ZEROES. 

INPUTS: 

VARO MATRIX TO BE CONDITIONED (N.N) 

IOPI PRINT OPTION; 0 NO PRINT, 1 PRINT 
N ACTUAL SIZE OF MATRIX VARO 

NMAX MAXIMUM SIZE OF N 

OUTPUTS : 

S CONDITIONED MATRIX (N.N) 

IOR BLOCK-DIAGONALIZING PERMUTATION INTEGER 

VECTOR (N) 

NBL INTEGER VECTOR OF SIZES OF EACH IRREDUCIBLE 
BLOCK (N) 

D VECTOR OF DIAGONAL ELEMENTS OF DIAGONAL SCALING 

MATRIX (N) 

TEMPORARY STORAGE: 


SS 

MATRIX (N.N) 

IN 

INTEGER 

VECTOR 

JBL 

INTEGER 

VECTOR 

IBL 

INTEGER 

VECTOR 

IC 

INTEGER 

VECTOR 


★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★A****************************** 


SUBROUTINE CONTRL (AA, BB, QC, NN, PCINV, KC, SS, CR, Cl, X, TS, 

1 XR, TT, AAA, EXT, AR, AI, IPER, IPERN, ADBLE, N, NC, N2, IOPI, 

2 IOP2, NMAX, NCMAX, N2MAX) 

* **★★★★***★★***★★*★★★*******************★★*★★***★*★★*★****★*★***★* 

SUBROUTINE CONTRL SOLVES THE OPTIMAL LINEAR REGULATOR PROBLEM. IT 
SETS UP AN N2 BY N2 MATRIX AAA, USING MATRICES AA, BB, QC, NN, AND 
PCINV. CONTRL OBTAINS THE SOLUTION TO THE RICCATI EQUATION, SS, 
AND THEN COMPUTES THE CONTROL GAINS, KC. CONTRL CALLS SUBROUTINES 
MATPRT AND RICSS. CONTRL IS CALLED BY SUBROUTINE AES800. 

INPUTS: 


AA SYSTEM MATRIX (N.N) 

BB CONTROL INPUT MATRIX (N.NC) 

QC STATE WEIGHTING MATRIX (N.N) 

NN STATE-CONTROL PRODUCT WEIGHTING MATRIX (N.NC) 

PCINV INVERSE OF CONTROL WEIGHTING MATRIX (NC.NC) 

IOPI SCALING PRINT OPTION: 0, NO PRINT; 1, PRINT 

I0P2 EIGENVECTOR PRINT OPTION: 0, NO PRINT; 1, PRINT 

N NUMBER OF STATE VARIABLES 

NC NUMBER OF CONTROL INPUTS 

N2 DIMENSION OF HAMILTONIAN MATRIX, 2 X N 

NMAX MAXIMUM SIZE OF N 

NCMAX MAXIMUM SIZE OF NC 

N2MAX MAXIMUM SIZE OF N2 



OUTPUTS: 


KC CONTROL GAIN MATRIX (NC.N) 

SS LQR RICCATI SOLUTION MATRIX (N,N) 

CR VECTOR OF REAL PARTS OF EIGENVALUES OF AAA (N2) 
Cl VECTOR OF IMAGINARY PARTS OF EIGENVALUES (N2) 

X MODIFIED EIGENVECTOR MATRIX OF AAA (N2,N2) 

TS SCALING TRANSFORMATION VECTOR OF AAA (N2) 

AAA HAMILTONIAN MATRIX FOR LQR RICCATI 
EQUATION (N2,N2) 

TEMPORARY STORAGE: 


XR 

MATRIX (N2,N2) 

TT 

MATRIX (N2,N2 

EXT 

MATRIX (N2,N2) 

AR 

VECTOR (N2) 

AI 

VECTOR (N2) 

IPER 

INTEGER VECTOR 

I PERN 

INTEGER VECTOR 

ADBLE 

VECTOR (N X N) 


★**★*★★★*****★**★**★*★**★***★*★**★★★★★*★**★*****★***★*★★★**★★**★*★ 


SUBROUTINE COVAR 
1 XX, YY, ZZ, UU, 



NM, NC 


, NO 




SUBROUTINE COVAR SETS UP MATRICES FOR SUBROUTINE LAPNV (LYAPUNOV 
EQUATION) WHICH IS THEN CALLED TO OBTAIN STATE COVARIANCE MATRIX, 
XX. XX, KALMAN FILTER ERROR COVARIANCE PP, AND CONTROL GAINS KC 
ARE USED TO OBTAIN CONTROL COVARIANCE - UU, OUTPUT COVARIANCE 
- YY, AND MEASUREMENT COVARIANCE - ZZ. COVAR CALLS SUBROUTINES 
ARRAY, LAPNV, AND MATPRT. COVAR IS CALLED BY SUBROUTINE AES800. 

INPUTS: 


AA SYSTEM MATRIX (N.N) 

BB CONTROL INPUT MATRIX (N,NC) 

HH MEASUREMENT MATRIX (NM,N) 

CC OUTPUT MATRIX (NO,N) 

DOUT FEED FORWARD MATRIX (NO.NC) 

QQ POWER SPECTRAL DENSITY MATRIX (N.N) 

(OF PLANT DISTURBANCE) 

PP KALMAN FILTER ERROR COVARIANCE MATRIX (N,N) 

KC CONTROL GAIN MATRIX (NC.N) 

N NUMBER OF STATE VARIABLES 

NM NUMBER OF MEASUREMENTS 

NC NUMBER OF CONTROL INPUTS 

NO NUMBER OF OUTPUTS 

NMAX MAXIMUM SIZE OF N 

NMMAX MAXIMUM SIZE OF NM 

NCMAX MAXIMUM SIZE OF NC 

NOMAX MAXIMUM SIZE OF NO 

OUTPUTS: 

XX STATE COVARIANCE MATRIX (N.N) 

YY OUTPUT COVARIANCE MATRIX (NO.NO) 

ZZ MEASUREMENT COVARIANCE MATRIX (NM,NM) 

UU CONTROL COVARIANCE MATRIX (NC.NC) 

TEMPORARY STORAGE: 

A MATRIX (N.N) 

Q MATRIX (N.N) 

WORK VECTOR (N) 


****************************************************************** 


SUBROUTINE CTBL (B, Cl, T, TINV, TINVB, EX1 , ADBLE, LEX, MEX, N 
1 NC, NMAX) 


SUBROUTINE CTBL COMPUTES THE (RELATIVE) CONTROLLABILITY OF A 
LINEAR SYSTEM DESCRIBED BY XDOT=A*X + B*U. 

NOTE: FOR A COMPLEX EIGENVALUE PAIR, THE CORRESPONDING TWO COLUMN 
ELEMENTS IN TINVB ARE STORED AS MAGNITUDE AND ANGLE (IN DEGREES) 
RESPECTIVELY. 

CTBL CALLS SUBROUTINES MATPRT AND MXINV. CTBL IS CALLED BY 
SUBROUTINE AES400. 

INPUTS: 

B SYSTEM INPUT MATRIX B (N,NC) 

Cl VECTOR OF IMAG PARTS OF THE EIGENVALUES (N) 

(OF MATRIX A) 

T MODIFIED EIGENVECTOR MATRIX OF MATRIX A (N,N) 

N NUMBER OF STATES 

NC NUMBER OF INPUTS 

NMAX MAXIMUM SIZE OF N 

OUTPUTS: 

TINV INVERSE OF MATRIX T (N,N) 

TINVB CONTROL EFFECTIVENESS MATRIX (N,NC) 

(IN MAGNITUDE AND PHASE ANGLE FORM) 

EX1 TINV*B WHERE TINV IS IN MODIFIED FORM (N,NC) 

TEMPORARY STORAGE: 

ADBLE VECTOR OF LENGTH N X N 

LEX INTEGER VECTOR (N) 

MEX INTEGER VECTOR (N) 

****************************************************************** 


SUBROUTINE DANSKY (A, X, Y, Z, N, NMAX) 
**★★★★**★★★★*★★★★★★***★★*******★★★★*****★★******★★★*★************* 

SUBROUTINE DANSKY COMPUTES THE COEFFICIENTS OF THE CHARACTERISTIC 
EQUATION. DANSKY CALLS SUBROUTINE POLMPY. DANSKY IS CALLED BY 
SUBROUTINE BOLLIN. 

INPUTS: 


AS CHARACTERISTIC EQUATION MATRIX (N,N) 

N ACTUAL SIZE OF MATRIX AS 

NMAX MAXIMUM SIZE OF N 

OUTPUTS: 

Z CHARACTERISTIC EQUATION COEFFICIENT VECTOR (N) 

TEMPORARY STORAGE: 

X VECTOR (N) 

Y VECTOR (N) 


***★**★★★★*★*★★★★★*★+**★********** ****************** ★★****★**+★*** 



SUBROUTINE DAVISO (A, B, C, N, NMAX, MC) 
****************************************************************** 

SUBROUTINE DAVISO TRANSFORMS X(DOT)=AX+BU, Y=C(TRANSPOSE)X USING 
Z=TX SUCH THAT Y IS A STATE VARIABLE OF Z(DOT)=TAT( I NVERSE ) +TBU - 
DAVISO DOES NOT CALL ANY SUBROUTINES. DAVISO IS CALLED BY 
SUBROUTINE BOLLIN. 


INPUTS: 


A 

SYSTEM MATRIX (N.N) 
SYSTEM VECTOR (N) 


B 


C 

SYSTEM VECTOR (N) 


N 

ACTUAL SIZE OF MATRIX A 


NMAX 

MAXIMUM SIZE OF N 

OUTPUTS: 


A 

TRANSFORMED MATRIX A (N, 


B 

TRANSFORMED VECTOR B (N) 


TEMPORARY STORAGE: 

MC INTEGER SCALAR 

****************************************************************** 


SUBROUTINE DSCA (A, R, CC, N, MS) 


SUBROUTINE DSCA FORMS R = A + CC * I, FOR EITHER VECTOR OR MATRIX 
IN VECTOR STORAGE MODE. DSCA DOES NOT CALL ANY SUBROUTINES. 

DSCA IS CALLED BY SUBROUTINE LAPNV. 


INPUTS: 

A INPUT MATRIX (N,N) , OR INPUT VECTOR (N) 

CC CONSTANT 

N ACTUAL SIZE OF SUBSCRIPT(S) OF MATRIX (VECTOR) A 

MS DECISION VARIABLE, 2 = MATRIX, 0 OR 1 = VECTOR 

OUTPUTS : 

R OUTPUT MATRIX (N.N), OR OUTPUT VECTOR (N) 


****************************************************************** 


SUBROUTINE DSCRT (DT, A, N, NMAX, ITIMES, EADT, INTGRL, C) 
****************************************************************** 

SUBROUTINE DSCRT CALCULATES EXP(A*DT) AND THE INTEGRAL 
FROM 0 TO DT OF EXP(A*T). AFTER EACH TERM OF THE SERIES IS 
MADE ON THE PERCENT CHANGE OCCURRING IN EACH TERM OF INTGRL. 

WHEN ALL CHANGES ARE LESS THAN .00001 , COMPUTATION IS STOPPED. 

IF ITIMES=50 BEFORE CONVERGENCE, DSCRT PROMPTS THE USER AS TO 
WHETHER TO COMPUTE MORE TERMS IN ORDER TO OBTAIN CONVERGENCE. 

DSCRT DOES NOT CALL ANY SUBROUTINES. DSCRT IS CALLED BY 
SUBROUTINE AES600. 


INPUTS 


DT TIME STEP 

A INPUT MATRIX (N,N) 

N ACTUAL SIZE OF MATRIX A 

NMAX MAXIMUM SIZE OF N 

OUTPUTS : 

ITIMES NO. OF TERMS IN SERIES EXPANSION 
EADT EXP(A*DT) (N,N) 

INTGRL INTEGRAL OF EXP(A*T) FROM T=0 TO T=DT (N,N) 
TEMPORARY STORAGE: 

C VECTOR (N) 


*★*★**★* ★★*★****************★*****************'**•**•*******•********* 


SUBROUTINE EGCK (AAA, X, CPR, CPI, EX1, EX2, PLAM, N, NMAX) 

****************************************************************** 

SUBROUTINE EGCK PERFORMS THE EIGENVALUE AND EIGENVECTOR CHECK. IT 
FORMS (AAA * X) AND (x * LAMBDA). EGCK DOES NOT CALL ANY 
SUBROUTINE. EGCK IS CALLED BY SUBROUTINE RICSS. 

INPUTS: 

AAA ORIGINAL MATRIX FOR WHICH EIGENVALUES 
AND EIGENVECTORS WERE FOUND (N,N) 

X MODIFIED EIGENVECTOR MATRIX (N,N) 

CPR VECTOR OF REAL EIGENVALUES (N) 

CPI VECTOR OF IMAGINARY EIGENVALUES (N) 

N ACTUAL SIZE OF MATRIX AAA 

NMAX MAXIMUM SIZE OF N 

OUTPUTS : 

EX 1 AAA * X MATRIX (N,N) 

EX2 X * LAMBDA MATRIX (N,N ) 

TEMPORARY STORAGE: 

PLAM MATRIX (N,N) 

****************************************************************** 


SUBROUTINE EGVCTR (AAA, CPR, CPI, X, N2, TT, EXT, AR, AI, IPERN, 

1 IPER, IOP2, N2MAX, ISEL, I HALF) 

**★*■*******★★★★★***★★★★**★★***★*****★■***************************** 

SUBROUTINE EGVCTR OBTAINS THE N2 BY N2 MODIFIEO EIGENVECTOR MATRIX 
X OF MATRIX AAA USING THE INVERSE ITERATION ALGORITHM. (THE 
EIGENVALUES OF AAA SHOULD HAVE BEEN PREVIOUSLY GENERATED USING 
SUBROUTINE EIGQR AND STORED IN CPR AND CPI.) IEND SPECIFIES THE 
NUMBER OF PASSES THRU THE INVERSE ITERATION ALGORITHM. 

IF ISEL=0, ALL EIGENVECTORS ARE OBTAINED. 

IF ISEL | 0 , ONLY THE ISEL (AND NEXT ONE IF A COMPLEX PAIR) IS 
OBTAINED. 

IF ISEL 0, THE ISELTH VECTOR IS PRINTED OUT AFTER EACH ITER. 

IF IHALF=1, THE FIRST N2/2 VECTORS ARE OBTAINED. 

IF I HALF .NE. 1, ALL VECTORS ARE OBTAINED. 

EGVCTR CALLS SUBROUTINES ARRAY, FACTR, MATPRT, AND PRMUTE. EGVCTR 
IS CALLED BY SUBROUTINES AES400, AES800, AND RICSS. 


INPUTS: 


AAA MATRIX FOR WHICH EIGENVECTORS ARE TO BE OBTAINED 
(N2,N2) 

CPR VECTOR OF REAL PARTS OF EIGENVALUES (N2) 

(OF AAA) 

CPI VECTOR OF IMAGINARY PARTS OF EIGENVALUES (N2) 

(OF AAA) 

N2 ACTUAL SIZE OF MATRIX AAA 

IOP2 PRINT OPTION: 0, NO PRINT; 1, PRINT 
N2MAX MAXIMUM SIZE OF N2 

ISEL SELECTION OPTION 

I HALF SELECTION OPTION 

OUTPUTS : 

X MODIFIED EIGENVECTOR MATRIX OF AAA (N2,N2) 

TEMPORARY STORAGE: 


TT 

MATRIX 

(N2,N2) 


EXT 

MATRIX 

(N2,N2) 


AR 

VECTOR 

(N2) 


AI 

VECTOR 

(N2) 


I PERN 

INTEGER 

VECTOR 

(N2) 

IPER 

INTEGER 

VECTOR (N2) 


****************************************************************** 


SUBROUTINE EIGEN (A, EIGR, EIGI, EX1, SSS, S, IA, IB, LEX, MEX, 

1 IBL, IC, EX4, N, NMAX) 

****************************************************************** 

SUBROUTINE EIGEN OBTAINS THE EIGENVALUES (EIGR AND EIGI) OF N X N 
MATRIX A BY FIRST REDUCING (IF A IS REDUCIBLE), THEN SCALING, 
HESSENBURG TRANSFORMING AND FINALLY APPLYING THE ' QR ' ALGORITHM 
EIGEN CALLS SUBROUTINES SCALEA, CONDI, ARRAY, HSBG, AND EIGQR. 
EIGEN IS CALLED BY SUBROUTINES AES400 AND AES800. 

INPUTS: 

A MATRIX TO GET EIGENVALUES FOR (N,N) 

N ACTUAL SIZE OF MATRIX A 

NMAX MAXIMUM SIZE OF N 

OUTPUTS: 

EIGR VECTOR OF REAL PARTS OF EIGENVALUES (N) 

EIGI VECTOR OF IMAGINARY PARTS OF EIGENVALUES (N) 

S MATRIX A IN REDUCED AND SCALED FORM (N.N) 

LEX BLOCK-DIAGONALIZING PERMUTATION INTEGER 

VECTOR (N) 

MEX INTEGER VECTOR OF SIZES OF EACH IRREDUCIBLE 

BLOCK (N) 

EX4 VECTOR OF DIAGONAL ELEMENTS OF DIAGONAL SCALING 

MATRIX (N) 

TEMPORARY STORAGE: 

EX 1 MATRIX (N.N) 

SSS MATRIX (N.N) 

I A INTEGER VECTOR (N) 

IB INTEGER VECTOR (N 

IBL INTEGER VECTOR (N) 

IC INTEGER VECTOR (N) 


****************************************************************** 



SUBROUTINE EIGQR (XR, N2, CR, Cl, IOP, N2MAX) 
****************************************************************** 

SUBROUTINE EIGQR COMPUTES THE EIGENVALUES OF MATRIX XR USING THE 
QR ALGORITHM. THIS MATRIX MUST BE IN UPPER HESSENBERG FORM. 

THE MAXIMUM NUMBER OF QR ITERATIONS USED IN FINDING ANY ONE 
EIGQR DOES NOT CALL ANY SUBROUTINES. EIGQR IS CALLED BY 
SUBROUTINES EIGEN, RICSS, AND ZEROES. 

INPUTS: 


XR MATRIX (IN UPPER HESSENBERG FORM) FOR WHICH 
EIGENVALUES ARE TO BE FOUND (N2,N2) 

N2 ACTUAL SIZE OF MATRIX XR 

IOP PRINT OPTION 

IOP 10, THE EIGENVALUES ARE WRITTEN ON UNIT 06 
I0P=0, NO WRITING TAKES PLACE 

IOP 0, THE EIGENVALUES ARE WRITTEN ON UNIT 06 AND 
ON UNIT 02 (TERMINAL) 

N2MAX MAXIMUM SIZE OF N2 

OUTPUTS: 

CR VECTOR OF REAL PARTS OF EIGENVALUES (N2) 

Cl VECTOR OF IMAGINARY PARTS OF EIGENVALUES (N2) 

****************************************************************** 


SUBROUTINE ESTMAT (AA, HH. QQ, RRINV, KE, PP, CR, Cl, X, TS, XR , 

1 TT, AAA, EXT, AR, AI, IPER, IPERN, ADBLE, N, NM, N2, IOPI, I0P2, 

2 NMAX, NMMAX, N2MAX) 

****************************************************************** 

SUBROUTINE ESTMAT SOLVES THE OPTIMAL LINEAR STATE ESTIMATION 
PROBLEM. IT SETS UP AN N2 BY N2 MATRIX AAA, USING MATRICES AA, 

HH, QQ, AND RRINV. ESTMAT OBTAINS THE KALMAN FILTER ERROR 
COVARIANCE, PP, AND THEN COMPUTES THE KALMAN FILTER GAINS, KE. 
ESTMAT CALLS SUBROUTINES MATPRT AND RICSS. ESTMAT IS CALLED BY 
SUBROUTINE AES800. 

INPUTS: 

AA SYSTEM MATRIX (N,N) 

HH MEASUREMENT MATRIX (NM,N ) 

QQ POWER SPECTRAL DENSITY MATRIX (N,N) 

(OF PLANT DISTURBANCE) 

RRINV INVERSE OF POWER SPECTRAL DENSITY MATRIX (NM,NM) 
(OF MEASUREMENT NOISE) 

IOPI SCALING PRINT OPTION: 0, NO PRINT; 1, PRINT 

I0P2 EIGENVECTOR PRINT OPTION: 0, NO PRINT; I, PRINT 

N NUMBER OF STATE VARIABLES 

NM NUMBER OF MEASUREMENTS 

N2 DIMENSION OF HAMILTONIAN MATRIX, 2 X N 

NMAX MAXIMUM SIZE OF N 

NMMAX MAXIMUM SIZE OF NM 

N2MAX MAXIMUM SIZE OF N2 

OUTPUTS : 

KE KALMAN FILTER GAIN MATRIX (N,NM) 

PP KALMAN FILTER ERROR COVARIANCE MATRIX (N,N) 


45 


CR 

VECTOR OF REAL PARTS OF EIGENVALUES (N2) 
(OF AAA) 

Cl 

VECTOR OF IMAGINARY PARTS OF EIGENVALUES (N2) 
(OF AAA) 

MODIFIED EIGENVECTOR MATRIX OF AAA (N2.N2) 

X 

TS 

SCALING TRANSFORMATION VECTOR OF AAA (N2) 

AAA 

HAMILTONIAN MATRIX FOR KALMAN FILTER RICCATI 
EQUATION (N2,N2) 

TEMPORARY STORAGE: 

XR 

MATRIX (N2.N2) 

TT 

MATRIX (N2,N2 

EXT 

MATRIX (N2,N2) 

AR 

VECTOR (N2) 

AI 

VECTOR (N2) 

IPER 

INTEGER VECTOR (N2) 

IPERN 

INTEGER VECTOR (N2) 

ADBLE 

VECTOR (N X N) 


* x-k'k'k'k'k'k'k-k'k'k+'k-k'k-k-k-k'k'k'k-k'k'k'k'irk-k'k'k'k'k'k'k'k'kirk-k-k-kick-k-k'kic'k'kicklrirk'k'ick'k'k'k'k’k'k'k-k 


SUBROUTINE FACTR (A, PER, N, IA, IER) 

★ ★***★★★**★★**★★★★★***★★★★★★**★*****★*★**★★★***★**★★★**★★★★★■*■*★*** 

SUBROUTINE FACTR FORMS THE LOWER AND UPPER TRIANGULAR MATRICES OF 
INPUT MATRIX A, SUCH THAT UPPER * LOWER = A. 

FACTR DOES NOT CALL ANY SUBROUTINES. FACTR IS CALLED BY 
SUBROUTINE EGVCTR. 

INPUTS: 

A INPUT MATRIX (N,N) 

N ACTUAL SIZE OF MATRIX A 

I A SAME AS N 

OUTPUTS: 

A INPUT MATRIX IN UPPER AND LOWER TRIANGULAR 

FORM (N,N) 

PER TRANSPOSITION VECTOR FOR MATRIX A (N) 

IER ERROR OPTION, 

IF IER ,NE. 0, FACTR IS WRONG 

IF IER .EQ. 0, FACTR HAS WORKED CORRECTLY 


SUBROUTINE FRPOLY (Zl, Z2, DD, HZ, G, AMP, PHA, N) 
*★★*★★★★★*★★*★★★*★***★*★★***★★■**★★★★★*★ *★**■**★★*★★★★*★★★★**★★★★★*★ 

SUBROUTINE FRPOLY EVALUATES TRANSFER FUNCTION Z2(S) / Z1(S) FOR 
S = 6.28 * HZ * J. FRPOLY DOES NOT CALL ANY SUBROUTINES. FRPOLY 
IS CALLED BY SUBROUTINE FRQP. 


INPUTS: 


Zl 

DENOMINATOR COEFFICIENT VECTOR 


12 

NUMERATOR COEFFICIENT VECTOR 


DD 

DOUT OR 0.0 


HZ 

FREQUENCY 


N 

SIZE OF COEFFICIENT VECTORS 

OUTPUTS: 


G 

COMPLEX TRANSFER FUNCTION 



AMP 

PHA 


AMPLITUDE 

PHASE 


****************************************************************** 


SUBROUTINE FRQP (Zl, Z2, DD, N, FI, DELF, IF, FREQ, AMP, PHASE, 

1 I SPACE, TSFTR) 

****************************************************************** 

SUBROUTINE FRQP GENERATES FREQUENCY RESPONSE AMP. AND PHASE, GIVEN 
TRANSFER FUNCTION NUMERATOR AND DENOMINATOR POLYNOMIAL 
COEFFICIENTS (GENERATED BY SUBROUTINE BOLLIN). 

FRQP CALLS SUBROUTINE FRPOLY, WHICH COMPUTES AMPLITUDE AND PHASE. 


FRQP IS 

CALLED BY SUBROUTINE FRSPNS. 

INPUTS: 

Zl 

DENOMINATOR POLYNOMIAL COEFFICIENT VECTOR 


ZZ 

NUMERATOR POLYNOMIAL COEFFICIENT VECTOR 


DD 

DOUT OR 0.0 


N 

SIZE OF COEFFICIENT VECTORS 


TSFTR 

TIME SCALE FACTOR 


FI 

INITIAL FREQUENCY 


DELF 

SPACING BETWEEN FREQUENCY POINTS 


IF 

NUMBER OF DESIRED POINTS TO BE GENERATED 


ISPACE 

CONTROLS FREQUENCY OF PRINTOUT FOR FREQ, 
AMP AND PHASE 

OUTPUTS 

FREQ 

FREQUENCY VECTOR 


AMP 

AMPLITUDE VECTOR 


PHASE 

PHASE VECTOR 


★★★★★★★★★★★★★★★★★★★★A********************************************* 


SUBROUTINE FRSPNS (A, B, C, DD, IOUT, JIN, N. TSFTR, DXM1, DXV1 , 

1 DXV2, DXV3, EXM1 , EXV1, EXV2, NMAX, NOMAX, DDCOF, DNCOF, FI, 

2 DELF, IF, FREQ, AMP, PHASE, ISPACE, IPRNT) 

****************************************************************** 

SUBROUTINE FRSPNS COMPUTES THE FREQUENCY RESPONSE OF THE IOUT 
OUTPUT TO THE JIN INPUT OF THE SYSTEM XDOT = A*X + B*U ; Y = C*X 
FRSPNS CALLS SUBROUTINES BOLLIN AND FRQP. FRSPNS IS CALLED BY 
SUBROUTINE AES500. 

INPUTS: 


A SYSTEM MATRIX (N,N) 

B SYSTEM MATRIX (N,NC) 

C SYSTEM MATRIX (NO,N) 

DD DOUT OR 0.0 

IOUT INDEX OF OUTPUT 

JIN INDEX OF INPUT 

N ACTUAL SIZE OF MATRIX A 

TSFTR TIME SCALE FACTOR 

NMAX MAXIMUM SIZE OF N 

NOMAX MAXIMUM SIZE OF NO 

FI INITIAL FREQUENCY 

DELF SPACING BETWEEN FREQUENCY POINTS 

IF NUMBER OF DESIRED POINTS TO BE GENERATED 

ISPACE CONTROLS FREQUENCY OF PRINTOUT FOR FREQ, 
AMP AND PHASE 

IPRNT PRINT OPTION, 0 IF STANDARD, 1 IF EXTENDED 



OUTPUTS: 


DXV3 

NUMERATOR COEFFICIENTS 

(N) 

EXM1 

A * TSFTR (N.N) 


EXV1 

JINTH ROW OF B * TSFTR 

(N) 

EXV2 

IOUTTH COLUMN OF C (N) 


DDCOF 

DENOMINATOR COEFFICIENTS (N) 

DNCOF 

NUMERATOR COEFFICIENTS (N) 

FREQ 

FREQUENCY VECTOR (500) 


AMP 

AMPLITUDE VECTOR (500) 


PHASE 

PHASE VECTOR (500) 



TEMPORARY STORAGE: 

DXM1 MATRIX (N,N) 

DXV1 VECTOR (N) 

DXV2 VECTOR (N) 

****************************************************************** 


SUBROUTINE GAIN (AA, BB, CC, DD, II, JJ, N, NZ, GAYN, EX1, EX4, - 
1 NMAX, NMMAX) 

****************************************************************** 

SUBROUTINE GAIN IS A COMPANION TO SUBROUTINE ZEROES. GAIN 
COMPUTES THE GAIN OF THE TRANSFER FUNCTION RELATING INPUT JJ AND 
OUTPUT II OF THE FOLLOWING NTH ORDER SYSTEM. 

(IN STATE VARIABLE FORM): 

XDOT = AA * X + BB * U 
Y=CC*X+DD*U 

GAIN DOES NOT CALL ANY SUBROUTINES. GAIN IS CALLED BY SUBROUTINE 
AES 700. 

INPUTS: 

AA SYSTEM MATRIX (N,N) 

BB CONTROL INPUT MATRIX (N, NUMBER OF POSSIBLE 

INPUTS) 

CC OUTPUT MATRIX (NUMBER OF POSSIBLE OUTPUTS, N) 

DD SCALAR RELATING U(JJ) TO Y( II ) 

II INDEX OF OUTPUT Y 

JJ INDEX OF INPUT U 

N ACTUAL NUMBER OF STATES 

NZ NUMBER OF NUMERATOR ZEROES IN TRANSFER FUNCTION 
(OBTAINED USING SUBROUTINE ZEROES) 

NMAX MAXIMUM SIZE OF N 

NMMAX MAXIMUM NUMBER OF OUTPUTS 

OUTPUTS: 

GAYN TRANSFER FUNCTION GAIN 

TEMPORARY STORAGE: 

EX1 MATRIX (N.N) 

EX4 VECTOR (N) 

****************************************************************** 


SUBROUTINE HSBG (N, A, IN) 

****************************************************************** 

SUBROUTINE HSBG REDUCES A MATRIX INTO UPPER ALMOST TRIANGULAR 
FORM. HSBG DOES NOT CALL ANY SUBROUTINES. HSBG IS CALLED BY 
SUBROUTINES EIGEN, RICSS, AND ZEROES. 
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INPUTS: 


N ACTUAL SIZE OF MATRIX A 

A INPUT MATRIX (N,N) 

IN MAXIMUM SIZE OF MATRIX A IN THE CALLING PROGRAM; 
IN = N, WHEN MATRIX A IS IN VECTOR STORAGE MODE. 

OUTPUTS : 

A OUTPUT MATRIX (N,N) 

****************************************************************** 

SUBROUTINE ICRSP (EX1, C, ICMTX, AMPIN, DT, TIME, TYOUT, XNEW, 

1 XOLD, TTIT, TTOP, TYTIT, IEXT, N, NOUT, NMAX, NOUTMX, ITRMX, IP, 

2 NAME, IONPLT) 

****************************************************************** 

SUBROUTINE ICRSP COMPUTES MULTIPLE INITIAL CONDITION RESPONSES OF 
THE SYSTEM: XDOT=A*X AND TYOUT=C*X 

BY SOLVING THE DIFFERENCE EQUATION: XNEW=EXl*XOLD 

THIS SUBROUTINE REQUIRES THAT THE STATE TRANSITION MATRIX, 
EXP(A*DT) , BE SUPPLIED AS INPUT MATRIX "EX1". DESIRED INITIAL 
CONDITION MAGNITUDES ARE SUPPLIED AS VECTOR 'AMPIN' AND THE 
DESIRED INITIAL CONDITION-OUTPUT RESPONSE COMBINATIONS ARE 
SELECTED BY APPROPRIATELY DEFINING ELEMENTS OF THE MATRIX 'ICMTX'. 
ICRSP CALLS PLOTTING SUBROUTINES ONLY. ICRSP IS CALLED BY 
SUBROUTINE AES600. 

INPUTS: 


EX 1 STATE TRANSITION, EXP(A*DT), MATRIX (N,N) 

C SYSTEM OUTPUT MATRIX (NOUT,N) 

ICMTX MATRIX OF ZEROES AND ONES (N,NOUT). 

ONES ARE PLACED IN SELECTED MATRIX POSITIONS TO 
INDICATE THE INITIAL CONDITION RESPONSES DESIRED. 
THE FIRST INDEX IS 'STATE', AND THE SECOND IS 
'OUTPUT'. THUS SUBROUTINE ICRSP MAY CALCULATE AS 
MANY AS N*NOUT INITIAL CONDITION RESPONSES. 

AMPIN VECTOR OF INPUT INITIAL CONDITION AMPLITUDES (N) 

DT TIME STEP 

TTIT PLOT TITLE (12) 

T TOP PLOT TITLE (12) 

TYTIT Y AXIS TITLE (4) 

N ACTUAL SIZE OF STATE TRANSITION MATRIX 

NOUT ACTUAL NUMBER OF POSSIBLE OUTPUTS 

NMAX MAXIMUM SIZE OF N 

NOUTMX MAXIMUM SIZE OF NOUT 

ITRMX NUMBER OF DESIRED TIME RESPONSE POINTS 
IP PLOT ENTITY INDEX (USED BY PLOTSUBS ONLY) 
INCREASES BY ONE FOR EACH FRAME 
NAME NAME OF PLOT DATASET (9) (USED BY PLOTSUBS ONLY) 
(PARTITIONED DATASET THAT HOLDS PLOT ENTITIES) 
IONPLT 0, IF OFFLINE PLOTS 
1, IF ONLINE PLOTS 

OUTPUTS: 

TIME VECTOR OF TIME POINTS (ITRMX) 

(SINGLE PRECISION) 

TYOUT MATRIX OF OUTPUT TRANSIENT RESPONSES FOR ANY 
ONE SPECIFIC INITIAL CGNTIDION. (ITRMX, NOUT) 
(SINGLE PRECISION) 

IP PLOT ENTITY INDEX (USED BY PLOTSUBS ONLY) 
INCREASES BY ONE FOR EACH FRAME 


TEMPORARY STORAGE: 

XNEW VECTOR (N) 

XOLD VECTOR (N) 

IEXT INTEGER VECTOR (N) 


SUBROUTINE LAPNV (A, X, B, QIN, NIN, WORK) 


****************************************************************** 

SUBROUTINE LAPNV SOLVES THE LYAPUNOV EQUATION, 

X*A‘ + A*X+B=0 

WHERE A 1 IS A TRANSPOSE 

A, B, AND X ARE ALL NXN ’MATRICES IN VECTOR STORAGE MODE, B IS 
SYMMETRIC ON INPUT, AND X IS SYMMETRIC ON OUTPUT. 

STEP 1 CALCULATE X(O) = A'*B*A 

STEP 2 THE EXACT SOLUTION X IS THE LIMIT OF THE SEQUENCE X(M) 
WHERE THE M REFERS TO THE M-TH TERM OF THE SEQUENCE. 

COMPUTE EACH TERM X(M+1) RECURSIVELY, BASED ON X(M) AS FOLLOWS, 
X(M+1 ) = X(M) + U(M)*X(M) * U'(M) 

U (0)= ( Q* I —A ' ) **(-1) * (Q*I+A‘ ) 

U(M) = U(O) **(2*M) 

LAPNV CALLS SUBROUTINES DSCA, MXINV, MXMLT, MXTRA, AND MXADD. 

LAPNV IS CALLED BY SUBROUTINES COVAR AND LYPCK. 

INPUTS: 

A LYAPUNOV EQUATION MATRIX (NIN, NIN) 

B LYAPUNOV EQUATION MATRIX (NIN, NIN) 

QIN CONVERGENCE FACTOR (TYPICALLY .1) 

NIN ACTUAL SIZE OF MATRIX A 

WORK(l) CONVERGENCE CHECK CRITERION (TYPICALLY l.E-6) 

OUTPUTS: 


X OUTPUT MATRIX (NIN, NIN) 

TEMPORARY STORAGE: 

WORK VECTOR (2 X NIN X NIN) 

****************************************************************** 


SUBROUTINE LOGSET (TMIN, TMAX, TNARR, INT) 
****************************************************************** 

SUBROUTINE LOGSET CALCULATES THE DECADES NECESSARY TO INCLUDE THE 
MINIMUM AND MAXIMUM OF THE DATA TO BE PLOTTED. 

LOGSET DOES NOT CALL ANY SUBROUTINES. 

INPUTS: 


TNARR ACTUAL DATA MINIMUM AND MAXIMUM (2) 
LOCATION 1 HOLDS THE MINIMUM 
LOCATION 2 HOLDS THE MAXIMUM 

OUTPUTS: 

TMIN MINIMUM DECADE FOR PLOTTING 

TMAX MAXIMUM DECADE FOR PLOTTING 

INT NUMBER OF INTERVALS FOR LOG AXIS 


****************************************************************** 


SUBROUTINE LYPCK (A, Q, XHT, E, R, EX1, WORK, N, NMAX) 


SUBROUTINE LYPCK COMPUTES THE RESIDUAL AND ERROR MATRICES 
ASSOCIATED WITH THE LYAPUNOV EQ. 

A*X + X*A**T + Q = 0 

WHERE 

SOLUTION— XHT 

RESIDUAL— R=A*XHT + XHT*A**T + Q 

ERROR E=XHT - X 

WHERE A*E + E*A**T - R = 0 

WORK VECTOR 'WORK 1 MUST BE DIMENSIONED |=2*N**2 
BEFORE COMPUTING E, R IS SYMMETRIZED, 

THE TRACES OF XHT, R, & E ARE PRINTED OUT; ALSO THE NORMALIZED 
ERROR INDEX, TR(E)/TR(XHT) , AND THE NORMALIZED DIAGONAL ELEMENTS 
OF THE ERROR MATRIX ARE PRINTED OUT. 

LYPCK CALLS SUBROUTINES MATPRT, ARRAY, AND LAPNV. LYPCK IS CALLED 
BY SUBROUTINE AES800. 

INPUTS: 

A LYAPUNOV EQUATION MATRIX (N,N) 

Q LYAPUNOV EQUATION MATRIX (N,N 

XHT LYAPUNOV SOLUTION MATRIX (N,N) 

N ACTUAL SIZE OF MATRIX A 

NMAX MAXIMUM SIZE OF N 

OUTPUTS: 

E ERROR MATRIX (N,N) 

R RESIDUAL MATRIX (N,N) 

TEMPORARY STORAGE: 

EX1 MATRIX (N,N) 

WORK VECTOR (2 X N X N) 

****************************************************************** 


SUBROUTINE MATCHG 

****************************************************************** 

SUBROUTINE MATCHG IS USED FOR CHANGING MATRICES AND DIMENSIONS 
USING NAMELIST ' MATDAT ' • THE CHANGES IN MATDAT ARE READ IN FROM 
THE TERMINAL. DATA IS TRANSFERRED TO PROGRAM AESOP VIA COMMONS 
■ABETC' AND 'DIMS 1 * MATCHG DOES NOT CALL ANY SUBROUTINES. MATCHG 
IS CALLED BY SUBROUTINE AES200. 

****************************************************************** 


SUBROUTINE MATIN 


SUBROUTINE MATIN IS USED FOR INPUTTING MATRICES AND DIMENSIONS 
FOR A 'SMALL* LQR /KALMAN FILTER PROBLEM TESTCASE. IT 
ALSO PRINTS OUT THE REFS NAMELIST. DATA IS PROVIDED FOR A 3RD 
ORDER TEST CASE HAVING TWO CONTROLS AND TWO SET-POINT OUTPUTS, TWO 
DISTURBANCES AND ONE NOISY MEASUREMENT. MATIN CALLS SUBROUTINE 
MATPRT. MATIN IS CALLED BY SUBROUTINE AES200. 



SUBROUTINE MATPRT (A, N, M, NMAX) 

****************************************************************** 

SUBROUTINE MATPRT PRINTS MATRIX A TEN COLUMNS PER PAGE, THE DEVICE 
ON WHICH THE PRINTING TAKES PLACE IS CONTROLLED BY 'IUNIT', 

I UN IT=2 — TERMINAL, IUNIT=6-~ LINEPRINTER. 

MATPRT DOES NOT CALL ANY SUBROUTINES. MATPRT IS CALLED BY 
SUBROUTINES AES200, AES300, AES400, AES600, AES700, AES800, COVAR, 
CTBL, EGCK, LYPCK, MATIN, MATRD , MODSHP, OBSBL, RESI, RICCHK, 
RICSS, AND UNRML. 

INPUTS: 


A MATRIX TO BE PRINTED (N,M) 

N NUMBER OF ROWS IN MATRIX A 

M NUMBER OF COLUMNS IN MATRIX A 

NMAX MAXIMUM SIZE OF N 

****************************************************************** 


SUBROUTINE MATRD 


SUBROUTINE MATRD IS USED FOR INPUTTING MATRICES, DIMENSIONS, AND 
REFS USING NAMELISTS 1 MATDAT 1 AND 'REFS'. MATDAT AND REFS ARE 
READ IN FROM UNIT 33. DATA IS TRANSFERRED TO PROGRAM AESOP VIA 
COMMONS ' COM1 1 , 'ABETC, 'DIMS', 'DIMS2', AND 'REF COM*. 

MATRD CALLS SUBROUTINE MATPRT. MATRD IS CALLED BY SUBROUTINE 
AES200. 

****************************************** ************************ 


SUBROUTINE MODSHP (A, B, Cl, N, NMAX) 
****************************************************************** 

SUBROUTINE MODSHP CALCULATES MODE SHAPES IN MAGNITUDE AND 
ANGLE (DEGREE) FORM. MODSHP CALLS SUBROUTINE MATPRT. MODSHP IS 
CALLED BY SUBROUTINES AES400, AES800, AND RICSS. 

INPUTS: 

A MODIFIED EIGENVECTOR MATRIX (N,N) 

Cl VECTOR OF IMAGINARY EIGENVALUES (N) 

N ACTUAL NUMBER OF EIGENVALUES 

NMAX MAXIMUM SIZE OF N 

OUTPUTS : 

B MODESHAPE IN MAGNITUDE AND ANGLE FORM (N,N) 

**************** ************************************************** 
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SUBROUTINE MXADD (A, B, R, N, M) 

****************************************************************** 


SUBROUTINE MXADD ADDS TWO IDENTICALLY SIZED MATRICES TO FORM A 
RESULTANT MATRIX. R CAN BE THE SAME AS A OR B IN THE CALLING 
PROGRAM. MXADD DOES NOT CALL ANY SUBROUTINES. MXADD IS CALLED 
BY SUBROUTINE LAPNV. 


INPUTS: 

A FIRST INPUT MATRIX (N,M) 

B SECOND INPUT MATRIX (N,M) 

N NUMBER OF ROWS IN MATRICES A, B, AND R 

M NUMBER OF COLUMNS IN MATRICES A, B, AND R 

OUTPUTS : 


R OUTPUT MATRIX (N # M) 


SUBROUTINE MXINV (A, N, D, L, M) 

* ************************* 

SUBROUTINE MXINV INVERTS A DOUBLE PRECISION MATRIX IN VECTOR 
STORAGE MODE BY USING THE GAUSS-JORDAN METHOD. THE DETERMINANT IS 
CALCULATED BUT IS ZERO IF THE MATRIX BEING INVERTED IS SINGULAR. 
MXINV DOES NOT CALL ANY SUBROUTINES. MXINV IS CALLED BY 
SUBROUTINES AES300, AES800, CTBL, LAPNV, AND RICSS. 

INPUTS: 

A MATRIX TO BE INVERTED, VECTOR STORAGE MODE (N,N) 

N ACTUAL SIZE OF MATRIX A 

OUTPUTS: 

A MATRIX INVERTED FORM, VECTOR STORAGE MODE (N,N) 

D SCALAR DETERMINANT (ZERO IF MATRIX A IS SINGULAR) 

TEMPORARY STORAGE: 

L INTEGER VECTOR (N) 

M INTEGER VECTOR (N) 

****************************************************************** 


SUBROUTINE MXMLT (A, B, R, N, L, M) 
****************************************************************** 

SUBROUTINE MXMLT MULTIPLIES TWO MATRICES IN VECTOR STOGAGE MODE 
TO FORM A RESULTANT MATRIX IN VECTOR STORAGE MODE. MXMLT DOES 
NOT CALL ANY SUBROUTINES. MXMLT IS CALLED BY SUBROUTINE LAPNV. 

INPUTS: 

A FIRST INPUT MATRIX (N,L) 

B FIRST INPUT MATRIX (L,M) 

N NUMBER OF ROWS IN MATRIX A 

L NUMBER OF COLUMNS IN MATRIX A 

AND ROWS IN MATRIX B 
M NUMBER OF COLUMNS IN MATRIX B 

OUTPUTS: 

R OUTPUT MATRIX (N,M) 



SUBROUTINE MXTRA (A, R, N, M) 


****************************************************************** 

SUBROUTINE MXTRA TRANSPOSES AN N BY M MATRIX A IN VECTOR STORAGE 
MODE TO FORM AN M BY N MATRIX R IN VECTOR STORAGE MODE. MXTRA 
DOES NOT CALL ANY SUBROUTINES. MXTRA IS CALLED BY SUBROUTINE 
LAPNV. 

INPUTS: 

A MATRIX TO BE TRANSPOSED (N,M) 

N NUMBERS OF ROWS IN MATRIX A 

AND COLUMNS IN MATRIX R 
M NUMBERS OF COLUMNS IN MATRIX A 

AND ROWS IN MATRIX R 

OUTPUTS: 

R RESULTANT MATRIX (M,N) 

★it**************************************************************** 


SUBROUTINE NRML (A, B, C, H, Q, RINV, D, DOUT, CSP, N, NC, NO, NM, 
1 NO, NMAX, NCMAX, NOMAX, NMMAX, FL34) 

****************************************************************** 


SUBROUTINE NRML READS FOUR NORMALIZATION VECTORS FROM NAMELIST 
NRMS AND NORMALIZES THE A, B, C, H, Q, RINV, D, DOUT, AND CSP 
MATRICES. THE SYSTEM THUS REPRESENTED IS DEFINED BY NORMALIZED 
STATE, CONTROL, OUTPUT, MEASUREMENT, AND SET POINT VECTORS. 

THE NORMALIZATION VECTORS ARE TRANSFERRED TO THE MAIN PROGRAM 
THROUGH COMMON 'NORMS'. NRML DOES NOT CALL ANY SUBROUTINES. 

NRML IS CALLED BY SUBROUTINE AES400. 

INPUTS: 

A UN-NORMALIZED SYSTEM MATRIX (N,N) 

B UN-NORMALIZED CONTROL INPUT MATRIX (N,NC) 

C UN-NORMALIZED OUTPUT MATRIX (NO,N) 

H UN-NORMALIZED MEASUREMENT MATRIX (NM,N) 

Q UN-NORMALIZED POWER SPECTRAL DENSITY MATRIX OF 

PLANT DISTURBANCE (N,N) 

RINV UN-NORMALIZED INVERSE OF POWER SPECTRAL DENSITY 
MATRIX OF MEASUREMENT NOISE (NM.NM) 

D UN-NORMALIZED DISTURBANCE INPUT MATRIX (N,ND) 

DOUT UN-NORMALIZED FEED FORWARD MATRIX FOR 
NON-ZERO SET POINT REGULATOR (NO.NC) 

CSP UN-NORMALIZED SET POINT OUTPUT MATRIX (NC,N) 

N ACTUAL NUMBER OF STATES 

NC ACTUAL NUMBER OF CONTROL INPUTS 

NO ACTUAL NUMBER OF OUTPUTS 

NM ACTUAL NUMBER OF MEASUREMENTS 

ND ACTUAL NUMBER OF DISTURBANCE INPUTS 

NMAX MAXIMUM SIZE OF N 

NCMAX MAXIMUM SIZE OF NC 

NOMAX MAXIMUM SIZE OF NO 

NMMAX MAXIMUM SIZE OF NM 

FL34 LOGICAL VARIABLE, ON INPUT 

TRUE, NORMALIZATION VECTOR INFORMATION 
(NAMELIST NRMS) HAS ALREADY BEEN READ IN 
FALSE, NORMALIZATION VECTOR INFORMATION 
(NAMELIST NRMS) NEEDS TO BE READ IN 
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OUTPUTS: 


A NORMALIZED SYSTEM MATRIX (N,N) 

B NORMALIZED CONTROL INPUT MATRIX (N,NC) 

C NORMALIZED OUTPUT MATRIX (NO,N) 

H NORMALIZED MEASUREMENT MATRIX (NM,N ) 

Q NORMALIZED POWER SPECTRAL DENSITY MATRIX OF 

PLANT DISTURBANCE (N,N) 

RINV NORMALIZED INVERSE OF POWER SPECTRAL DENSITY 
MATRIX OF MEASUREMENT NOISE (NM.NM) 

D NORMALIZED DISTURBANCE INPUT MATRIX (N,ND) 


DOUT NORMALIZED FEED FORWARD MATRIX FOR 
NON-ZERO SET POINT REGULATOR (NO.NC) 

CSP NORMALIZED SET POINT OUTPUT MATRIX (NC,N) 

FL34 LOGICAL VARIABLE, ON OUTPUT SET TO 

TRUE IF NAMELIST NRMS HAS BEEN READ IN 


SUBROUTINE OBSBL (H, T, Cl, HT, EX1, N, NR, NRMAX, NMAX) 


SUBROUTINE OBSBL COMPUTES THE OBSERVABILITY MATRIX HT FOR THE 
LINEAR SYSTEM DESCRIBED BY XDOT=A*X + B*U, AND Y=H*X. 

NOTE: FOR A COMPLEX EIGENVALUE PAIR, THE CORRESPONDING TWO COLUMN 
ELEMENTS IN HT ARE STORED AS MAGNITUDE AND ANGLE (IN DEGREES) 
RESPECTIVELY. 

OBSBL CALLS SUBROUTINE MATPRT. OBSBL IS CALLED BY SUBROUTINE 


AES400. 

INPUTS: 


H 

SYSTEM OUTPUT MATRIX (NR,N) 


T 

MODIFIED EIGENVECTOR MATRIX OF MATRIX A (N,N) 


Cl 

VECTOR OF IMAG PARTS OF THE EIGENVALUES (N) 


N 

(OF MATRIX A) 

ACTUAL NUMBER OF COLUMNS IN MATRIX H 


NR 

ACTUAL NUMBER OF ROWS IN MATRIX H 


NRMAX 

MAXIMUM SIZE OF NR 


NMAX 

MAXIMUM SIZE OF N 

OUTPUTS: 


HT 

OBSERVABILITY MATRIX (NR,N) 


EX1 

(IN MAGNITUDE AND PHASE ANGLE FORM) 
OBSERVABILITY MATRIX (NR,N) 


****************************************************************** 


SUBROUTINE ORDER (CR, Cl, NE, EPS) 


GIVEN A SET OF NE EIGENVALUES, SYMMETRICALLY LOCATED WITH RESPECT 
TO THE IMAGINARY AXIS, SUBROUTINE ORDER PLACES ONES WITH POSITIVE 
REAL PARTS IN FIRST NE/2 LOCATIONS. CORRESPONDING SYMMETRIC 
EIGENVALUES WITH NEGATIVE REAL PARTS ARE PUT IN LOCATIONS 
NE/2 + 1 THROUGH NE. EPS IS THE CONVERGENCE CRITERION USED IN 
DETERMINING IF A PAIR OF EIGENVALUES ARE SYMMETRIC. ORDER DOES 
NOT CALL ANY SUBROUTINES. ORDER IS CALLED BY SUBROUTINE RICSS. 
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INPUTS: 

CR VECTOR OF REAL PARTS OF EIGENVALUES 
UNORDERED (NE) 

Cl VECTOR OF IMAGINARY PARTS OF EIGENVALUES 
UNORDERED (NE) 

NE NUMBER OF EIGENVALUES 

EPS CRITERION FOR SYMMETRY 

OUTPUTS: 

CR VECTOR OF REAL PARTS OF EIGENVALUES 

ORDERED (NE) 

Cl VECTOR OF IMAGINARY PARTS OF EIGENVALUES 

ORDERED (NE) 

****************************************************************** 


SUBROUTINE POLMPY (X, Y, Z, NX, NY, NZ) 
******************************************************************* 

SUBROUTINE POLMPY MULTIPLIES X*Y=Z (THE LEADING POLYNOMIAL 
COEFFICIENT IS ASSUMED TO BE UNITY). POLMPY DOES NOT CALL ANY 
ANY SUBROUTINES. POLMPY IS CALLED BY SUBROUTINE DANSKY. 

INPUTS: 

X POLYNOMIAL COEFFIENT VECTOR (NX) 

Y POLYNOMIAL COEFFIENT VECTOR (NY) 

NX ORDER OF POLYNOMIAL FOR WHICH VECTOR X IS THE 
LIST OF COEFFICIENTS (OTHER THAN THE FIRST) 

NY ORDER OF POLYNOMIAL FOR WHICH VECTOR Y IS THE 
LIST OF COEFFICIENTS (OTHER THAN THE FIRST) 

OUTPUTS: 

Z X VECTOR * Y VECTOR (NZ) 

NZ ORDER OF POLYNOMIAL FOR WHICH VECTOR Z IS THE 

LIST OF COEFFICIENTS (OTHER THAN THE FIRST) 


SUBROUTINE PREREQ (WHEN, ICA, I AND, MZ, II) 


SUBROUTINE PREREQ CHECKS TO SEE IF PREREQUISITE CASES HAVE BEEN 
DONE FOR THE CASE ABOUT TO BE RUN. IF NOT, IT PRINTS OUT WHAT THE 
PREREQUISITES ARE. PREREQ DOES NOT CALL ANY SUBROUTINES. PREREQ 
IS CALLED BY SUBROUTINES AESIOO, AES200, AES300, AES400, AES500, 
AES600, AES 700, AND AES800. 

INPUTS: 

WHEN LOGICAL MATRIX OF PREREQS (450,50) 

ICA VECTOR OF CASE NUMBERS TO BE DONE (1000) 

MZ WHICH CASE IS TO BE CHECKED FOR PREREQUISITES 

II WHICH ROW OF MATRIX 'WHEN' IS TO BE LOOKED AT 

OUTPUTS: 

I AND DECISION VARIABLE, 0 IF PREREQUISITES HAVE BEEN 

DONE, 1 IF PREREQUISITES HAVE NOT BEEN DONE 

****************************************************************** 
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SUBROUTINE PRMUTE (X, ITRANS, IA, NE, LA, NEMAX) 


SUBROUTINE PRMUTE PERMUTES ELEMENTS IN LA COLUMN OF NE BY NE 
MATRIX X AS DICTATED BY TRANSPOSITION VECTOR ITRANS. ITRANS IS 
PRODUCED, IN THIS CASE, BY SUBROUTINE FACTR. 

PRMUTE DOES NOT CALL ANY SUBROUTINES. PRMUTE IS CALLED BY 
SUBROUTINE EGVCTR. 

INPUTS: 

X MATRIX TO BE PERMUTED ( NE, NE) 

ITRANS TRANSPOSITION VECTOR (NE) 

NE ACTUAL SIZE OF MATRIX X 

LA SPECIFIC COLUMN OF MATRIX X TO BE PERMUTED BY 
THE TRANSPOSITION VECTOR ITRANS 
NEMAX MAXIMUM SIZE OF NE 

OUTPUTS : 

X INPUT MATRIX IN PERMUTED FORM (NE.NE) 

TEMPORARY STORAGE: 

IA INTEGER VECTOR (NE) 

****************************************************************** 


SUBROUTINE REDU (VARO, SS, S, IN, JBL, INBL, IOR, NBL, IBL, IC, N, 
1 NMAX) 

****************************************************************** 

SUBROUTINE REDU USES HARARYS METHOD FOR REDUCTION OF A REDUCIBLE 
MATRIX TO BLOCK DIAGONAL FORM. REDU DOES NOT CALL ANY SUBROUTINE. 
REDU IS CALLED BY SUBROUTINE CONDI. 


INPUTS: 


VARO 

MATRIX TO BE REDUCED (N,N) 


N 

ACTUAL SIZE OF MATRIX VARO 


NMAX 

MAXIMUM SIZE OF N 


OUTPUTS: 

S 

BLOCK DIAGONAL MATRIX (N,N) 


INBL 

NUMBER OF REDUCIBLE BLOCKS 


IOR 

BLOCK-DIAGONALIZING PERMUTATION 
VECTOR (N) 

INTEGER 

NBL 

INTEGER VECTOR OF SIZES OF EACH 
BLOCK (N) 

IRREDUCIBLE 


TEMPORARY STORAGE: 


SS 

MATRIX (N,N) 

IN 

INTEGER 

VECTOR 

JBL 

INTEGER 

VECTOR 

IBL 

INTEGER 

VECTOR 

IC 

INTEGER 

VECTOR 


II 


II I III III I III INI I I III i i 
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SUBROUTINE RESI (C, B, EIGR, EIGI, R, EXA, N, NC, NO, NMAX, NOMAX) 
****************************************************************** 

SUBROUTINE RESI COMPUTES THE RESIDUE MATRICES FOR THE LINEAR 

SYSTEM, XDOT=A*X + B*U, AND, Y=C*X 

WHERE THE SYSTEM IS ASSUMED TO BE IN BLOCK-DIAGONAL FORM. 

MATRICES C AND B ARE INPUT TO THE PROGRAM. MATRIX C IS ASSUMED TO 
HAVE BEEN TRANSFORMED TO THE FORM CORRESPONDING TO BLOCK-DIAGONAL 
A MATRIX USING SUBROUTINE OBSBL. MATRIX B IS ASSUMED TO HAVE BEEN 
SIMILARLY TRANSFORMED USING SUBROUTINE CTBL. 

NOTE: TWO RESIDUE MATRICES ARE PRINTED OUT FOR A COMPLEX 

EIGENVALUE; THE FIRST CONTAINS RESIDUE MAGNITUDES AND THE SECOND 
CONTAINS RESIDUE PHASE ANGLES IN DEGREES. 

RESI CALLS SUBROUTINE MATPRT. RESI IS CALLED BY SUBROUTINE AES400. 
INPUTS: 

C OUTPUT MATRIX (NO,N) 

B INPUT MATRIX (N,NC) 

EIGR VECTOR OF REAL PARTS OF EIGENVALUES (N) 

EIGI VECTOR OF IMAGINARY PARTS OF EIGENVALUES (N) 

N NUMBER OF STATES 

NC NUMBER OF INPUTS 

NO NUMBER OF OUTPUTS 

NMAX MAXIMUM SIZE OF N 

NOMAX MAXIMUM SIZE OF NO 

OUTPUTS: 

R RESIDUE ARRAY (N,NO,NC) 

TEMPORARY STORAGE: 

EXA MATRIX (NO.NC) TEMPORARILY STORES EACH RESIDUE 

MATRIX BEFORE BEING PRINTED OUT 


SUBROUTINE RICCHK (AAA, S, R, N, NMAX, N2MAX ) 
****************************************************************** 

SUBROUTINE RICCHK COMPUTES THE RESIDUAL ERROR MATRIX FOR THE 
RICCATI EQUATION 

S*AAA(22)*T + AAA(22)*S - S*AAA(12)*S + AAA(21) = 0 
OR, -S*AAA(11) - AAA(11 )*T*S - $*AAA(12)*S + AAA(21) = 0 

WHERE 

S =RICCATI SOLUTION MATRIX 

AAA=THE FOUR N BY N BLOCKS OF THE HAMILTONIAN MATRIX 
R “RES I DUAL ERROR MATRIX 
R IS GIVEN BY 

R“ S*AAA(22)*T + AAA(22)*S - S*AAA(12)*S + AAA(21) 

=-(S*AAA(n) + AAA(11)*T*S - S*AAA(12)*S - AAA(21J) 

RICCHK CALLS SUBROUTINE MATPRT. RICCHK IS CALLED BY SUBROUTINE 
AES800. 

INPUTS: 


AAA HAMILTONIAN MATRIX (2 X N,2 X N) 

S RICCATI SOLUTION MATRIX (N,N) 

N ACTUAL SIZE OF MATRIX S 

NMAX MAXIMUM SIZE OF N 

N2MAX 2 X NMAX 
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OUTPUTS: 


R RESIDUAL ERROR MATRIX (N,N) 

****************************************************************** 


SUBROUTINE RICSS (AAA, X, OUTPUT, CR, Cl, TS, XR, EXT, TT, IPER, 

1 IPERN, AR, AI, ADBLE, IOP1, I0P2, N, N2, NMAX, N2MAX) 

****************************************************************** 

SUBROUTINE RICSS COMPUTES THE OUTPUT SOLUTION TO THE STEADY STATE 
MATRIX RICCATI EQUATION. THE INPUT IS AN N2 BY N2 MATRIX, AAA, 
WHICH IS THE HAMILTONIAN MATRIX FOR KALMAN FILTER MATRIX RICCATI 
EQUATION. RICSS CALLS SUBROUTINES ARRAY, MXINV, EGCK, EGVCTR, 
EIGQR, HSBG, MATPRT, MODSHP, ORDER, AND SCALEA. RICSS IS CALLED 
BY SUBROUTINES CONTRL AND ESTMAT. 

INPUTS: 

AAA HAMILTONIAN MATRIX FOR KALMAN FILTER RICCATI 
EQUATION (N2,N2) 

IOP1 SCALING PRINT OPTION: 0, NO PRINT; 1, PRINT 

I0P2 EIGENVECTOR PRINT OPTION: 0, NO PRINT; 1, PRINT 

N NUMBER OF STATE VARIABLES 

N2 DIMENSION OF HAMILTONIAN MATRIX, 2 X N 

NMAX MAXIMUM SIZE OF N 

N2MAX MAXIMUM SIZE OF N2 

OUTPUTS : 

X MODIFIED EIGENVECTOR MATRIX OF AAA (N2.N2) 

OUTPUT RICCATI SOLUTION MATRIX (N,N) 

CR VECTOR OF REAL PARTS OF EIGENVALUES (N2) 

(OF AAA) 

Cl VECTOR OF IMAGINARY PARTS OF EIGENVALUES (N2) 

(OF AAA) 

TS SCALING TRANSFORMATION VECTOR OF AAA (N2) 
TEMPORARY STORAGE: 


XR 

MATRIX l 

(N2,N2) 

EXT 

MATRIX ( 

:N2,N2 

TT 

MATRIX i 

[N2,N2) 

IPER 

INTEGER 

VECTOR 

IPERN 

INTEGER 

VECTOR 

AR 

VECTOR | 

:N2) 

AI 

VECTOR < 

N2) 

ADBLE 

VECTOR ( 

;n x n) 


•kit'k-kirk'k'klrk'kicirkif'k-k-kitir'k'k'kirk'k'kirk'kie'kirk-kic'k'kirkicie'kirirkicic'k-kititit'k-k'k'kicirk'kititit'k'k 


SUBROUTINE SCALEA (A, TS, N2, I0P1, N2MAX) 
****************************************************************** 

SUBROUTINE SCALEA TRANSFORMS N2 BY N2 MATRIX A USING DIAGONAL 
MATRIX TS SO THAT THE NORM OF A IS MINIMIZED. THE RESULTING SCALED 
MATRIX IS STORED IN A. IF SCALEA FINDS A TO BE REDUCIBLE, IER IS 
SET TO 1. SCALEA DOES NOT CALL ANY SUBROUTINES. SCALET IS CALLED 
BY SUBROUTINES CONDI, EIGEN, AND RICSS. 

INPUTS: 


A MATRIX TO BE SCALED (N2.N2) 

N2 ACTUAL SIZE OF MATRIX A 

I0P1 PRINT OPTION; 0 NO PRINT, 1 PRINT 

N2MAX MAXIMUM SIZE OF N2 


OUTPUTS: 

A INPUT MATRIX IN SCALEO FORM (N2,N2) 

TS VECTOR OF DIAGONAL ELEMENTS OF DIAGONAL SCALING 
MATRIX (N2) 

* ***★********★***•*★*****★***★***★*★******★***★★******★********★★** 


SUBROUTINE STP (EX1, EX2, B, C, DOUT, IOMTX, AMPIN, DT, TIME, 

1 TYOUT, XNEW, XOLD, ANR, TTIT, TTOP, TYTIT, IEXT, N, NIN, NOUT, 

2 NMAX, NINMAX, NOUTMX, ITRMX, IP, NAME, IONPLT) 

****************************************************************** 

SUBROUTINE STP COMPUTES MULTIPLE STEP RESPONSES OF THE SYSTEM 
XDOT=A*X+B*U; TYOUT=C*X+DOUT*U 
BY SOLVING THE DIFFERENCE EQ. 

XNEW=EXl*XOLD + EX2*B*AMPIN(L) 

THIS SUBROUTINE REQUIRES THAT THE STATE TRANSITION MATRIX, 
EXP(A*DT) , AND ITS INTEGRAL FROM TIME=0 TO TIME=DT, BE SUPPLIED AS 
INPUT MATRICES 'EXT AND 'EX2' . DESIRED INPUT STEP MAGNITUDES ARE 
SUPPLIED AS VECTOR 'AMPIN' AND THE DESIRED STEP INPUT-OUTPUT 
RESPONSE COMBINATIONS ARE SELECTED BY APPROPRIATELY DEFINING 
ELEMENTS OF THE MATRIX 'IOMTX'. STP CALLS PLOTTING SUBROUTINES 
ONLY. STP IS CALLED BY SUBROUTINE AES600. 

INPUTS: 


EX1 STATE TRANSITION, EXP(A*DT), MATRIX (N,N) 

EX2 INTEGRAL OF THE STATE TRANSITION MATRIX FROM 
TIME=0 TO TIME=DT (N,N) 

B (CONTINUOUS) SYSTEM INPUT MATRIX (N,NIN) 

C SYSTEM OUTPUT MATRIX (NOUT.N) 

DOUT SYSTEM INPUT/OUTPUT FEEDTHRU MATRIX (NOUT.N IN) 
IOMTX MATRIX OF ZEROES AND ONES (NIN, NOUT). 

ONES ARE PLACED IN SELECTED MATRIX POSITIONS TO 
INDICATE THE STEP RESPONSES DESIRED. THE FIRST 
INDEX IS 'INPUT', THE SECOND IS ’OUTPUT'. THUS 
SUBROUTINE STP MAY CALCULATE AS MANY AS NIN*NOUT 
STEP RESPONSES. 

AMPIN VECTOR OF INPUT STEP AMPLITUDES (NIN) 

DT TIME STEP 
TTIT PLOT TITLE (12) 

TTOP PLOT TITLE (12) 

TYTIT Y AXIS TITLE (4) 

N ACTUAL SIZE OF STATE TRANSITION MATRIX 

NIN ACTUAL NUMBER OF POSSIBLE INPUTS 

NOUT ACTUAL NUMBER OF POSSIBLE OUTPUTS 

NMAX MAXIMUM SIZE OF N 
NINMAX MAXIMUM SIZE OF NIN 
NOUTMX MAXIMUM SIZE OF NOUT 
ITRMX NUMBER OF DESIRED TIME RESPONSE POINTS 
IP PLOT ENTITY INDEX (USED BY PLOTSUBS ONLY) 
INCREASES BY ONE FOR EACH FRAME 
NAME NAME OF PLOT DATASET (9) (USED BY PLOTSUBS ONLY) 
(PARTITIONED DATASET THAT HOLDS PLOT ENTITIES) 
IONPLT 0, IF OFFLINE PLOTS 
1, IF ONLINE PLOTS 

OUTPUTS: 

TIME VECTOR OF TIME POINTS (ITRMX) 

(SINGLE PRECISION) 

TYOUT MATRIX OF OUTPUT TRANSIENT RESPONSES FOR ANY 
ONE SPECIFIC INPUT STEP (ITRMX, NOUT) 

(SINGLE PRECISION) 

IP PLOT ENTITY INDEX (USED BY PLOTSUBS ONLY) 
INCREASES BY ONE FOR EACH FRAME 



TEMPORARY STORAGE 


XNEW VECTOR (N) 

XOLD VECTOR ( N ) 

ANR VECTOR (N) 

IEXT INTEGER VECTOR (N) 


****************************************************************** 


SUBROUTINE UNRML (KC, KE, KFF, PP, NC, NM, N, NCMAX, NMAX, FL34) 
****************************************************************** 

SUBROUTINE UNRML IS USED TO CONVERT NORMALIZED KC, KE, KFF, & PP 
MATRICES TO UN-NORMALIZED FORM. NORMALIZATION VECTOR INFORMATION 
IS FED IN THRU COMMON 'NORMS* OR IF NECESSARY, READ IN OFF UNIT 34 
AS NAMELIST NRMS. IF FL34 IS .TRUE., IT MEANS THAT NRMS HAS 
ALREADY BEEN READ IN BY SUBROUTINE NRML & THUS IT SHOULDN'T BE 
READ IN HERE. UNRML CALLS SUBROUTINE MATPRT. UNRML IS CALLED BY 
SUBROUTINE AES400. 

INPUTS: 

KC NORMALIZED CONTROL GAIN MATRIX (NC,N) 

KE NORMALIZED KALMAN FILTER GAIN MATRIX (N,NM) 

KFF NORMALIZED FEED FORWARD GAIN MATRIX FOR 
NON-ZERO SET POINT REGULATOR (NC,NC) 

PP NORMALIZED KALMAN FILTER ERROR COVARIANCE 
MATRIX (N,N) 

NC ACTUAL NUMBER OF CONTROL INPUTS 

NM ACTUAL NUMBER OF MEASUREMENTS 

N ACTUAL NUMBER OF STATES 

NCMAX MAXIMUM SIZE OF NC 

NMAX MAXIMUM SIZE OF N 

FL34 LOGICAL VARIABLE, ON INPUT 

TRUE, NORMALIZATION VECTOR INFORMATION 
{NAMELIST NRMS) HAS ALREADY BEEN READ IN 
FALSE, NORMALIZATION VECTOR INFORMATION 
(NAMELIST NRMS) NEEDS TO BE READ IN 

OUTPUTS : 

KC UN-NORMALIZED CONTROL GAIN MATRIX (NC,N) 

KE UN-NORMALIZED KALMAN FILTER GAIN MATRIX (N,NM) 

KFF UN-NORMALIZED FEED FORWARD GAIN MATRIX FOR 
NON-ZERO SET POINT REGULATOR (NC,NC) 

PP UN-NORMALIZED KALMAN FILTER ERROR COVARIANCE 

MATRIX (N,N) 

FL34 LOGICAL VARIABLE, ON OUTPUT SET TO 

TRUE IF NORMALIZATION VECTOR INFORMATION 
(NAMELIST NRMS) HAS BEEN READ IN 

****************************************************************** 


SUBROUTINE UZR901 

****************************************************************** 
THIS IS FOR A USER-WRITTEN ROUTINE 


****************************************************************** 



SUBROUTINE UZR902 

****************************************************************** 

THIS IS FOR A USER-WRITTEN ROUTINE 
****************************************************************** 


SUBROUTINE UZR903 

****************************************************************** 

THIS IS FOR A USER-WRITTEN ROUTINE 
****************************************************************** 


SUBROUTINE UZR904 

****************************************************************** 

THIS IS FOR A USER-WRITTEN ROUTINE 
****************************************************************** 


SUBROUTINE ZEROES (AA, BB, CC, DD, CONST, ANR, ANI, N, II, JJ, L, 

1 AR, TS, BR, LWV, MWV, ZERMAX, IA, IB, IBL, IC, S, SS, NMAX, 

2 NOMAX) 

****************************************************************** 

SUBROUTINE ZEROES FINDS THE NUMERATOR ZEROES OF THE TRANSFER 
FUNCTION Y( II) / U(JJ) = CC * ( (S * I - A ) ** -1 ) * BB 
II DENOTES DESIRED COMPONENT OF OUTPUT VECTOR 
JJ DENOTES DESIRED COMPONENT OF INPUT VECTOR 

ZEROES CALLS SUBROUTINES CONDI, EIGQR, HSBG, AND ARRAY. ZEROES IS 
CALLED BY SUBROUTINE AES700. 

INPUTS: 

AA SYSTEM MATRIX (N,N) 

BB INPUT MATRIX (N, NUMBER OF POSSIBLE INPUTS) 

CC OUTPUT MATRIX (NUMBER OF POSSIBLE OUTPUTS, N) 

DD CONSTANT 

CONST ITERATION CONSTANT 

N ACTUAL SIZE OF MATRIX AA 

II OUTPUT COMPONENT 

JJ INPUT COMPONENT 

NMAX MAXIMUM SIZE OF N 

NOMAX MAXIMUM NUMBER OF OUTPUTS 

OUTPUTS: 

ANR VECTOR OF REAL PARTS OF EIGENVALUES (N) 

ANI VECTOR OF IMAGINARY PARTS OF EIGENVALUES (N) 

L NUMBER OF ZEROES 

ZERMAX MAXIMUM EXPECTED VALUE OF TRANSFER FUNCTION 
ZEROES 



TEMPORARY STORAGE: 


AR 

TS 

BR 

MATRIX 

MATRIX 

VECTOR 

[N,N) 

Sr 


LWV 

INTEGER 

VECTOR ( 

.2 X N) 

MWV 

INTEGER 

VECTOR 

i2 X N) 

IA 

INTEGER 

VECTOR | 

2 X N 

IB 

INTEGER 

VECTOR i 

(2 X N) 

IBL 

INTEGER 

VECTOR | 

'2 X N 

IC 

S 

SS 

INTEGER 
MATRIX | 
MATRIX 1 

VECTOR i 

;n,n 

(N,N) 

[2 X N) 


* ★★★★★*★★★***★★*★*★★***★ ★★★*★*★*★**★★★★★*****************★** ****** 



Appendix C 
Test Cases 


This appendix presents the results of two test cases. Test case I exercises almost all of the 77 
available AESOP functions. Test case II shows how AESOP can be used in an interactive manner to 
design a control system. 

Test Case I - Third-Order Problem to Demonstrate Full AESOP Program Capabilities 

A block diagram of the selected third-order, open-loop plant is shown in figure 14. The plant is 
characterized by three states, two controls, one noisy measurement, two plant-noise disturbances, 
two outputs, and two set -point variables. The plant is stable with one real pole at 0.1 rad/sec and a 
complex pole pair with a natural frequency of 1.0 rad/sec and a damping ratio of 0.001. The com- 
putations performed do not all relate to one specific design problem but are done so as to exercise all 
of the AESOP functions. For example, frequency responses are computed for a system with perfect 
state measurement (x!,x 2 ,X 3 ) and state feedback as well as for Kalman filter feedback using a single 
noisy measurement, z\. 

The input matrices for test case I are generated by function 201 and are listed in table V. They 
include the plant matrices A, B, C, D, CSP, H, and DOUT as well as quadratic weighting matrices 
QC, NN, and PCINV and noise matrices QQ and RRINV. The test case exercises the AESOP 
functions shown in table XI. 

Pages 65 to 80 are the terminal printout sheets produced at the user’s terminal. Notes have been 
added to indicate what results are being generated and displayed. In the pocket at the rear of this 
report are microfiche copies of the output dataset and plots generated for this test case. The output 
dataset contains the complete results generated by the AESOP run, including results that were 
printed out at the user’s terminal. The CPU time required to run this test case on an IBM 370/3033 is 
approximately 30 to 40 seconds depending on which device was used to produce graphic output. 



Plant noise vector 




measurement noise » 


Output vector ^y^); set-point vector “(y^) 


Figure 14. - Block diagram of third-order system used for test case I. 





TABLE XI. - TASKS PERFORMED IN TEST CASE I 


Design task performed 

Page number, 
appendix C 

Open-loop analyses 

68 

Optimal regulator gains and associated error checks 

68 

Kalman filter gains and associated error checks 

70 

State covariance matrices and feedforward matrix 

72 

Frequency responses (and Bode plots) 

72 

Transient responses (and plots) 

74 

Transfer function gain and zeros 

76 

System matrix normalization 

77 

Repeat of tasks 2 and 3 with normalized matrices 

79 

Unnormalization of resultant matrices 

79 


Terminal Printout for Test Case I 


AES RUN 1 User invokes PROCDEF AESRUN to begin run 

CANCELLED: 0UT1 UNKNOWN. 

0UT1 IS DDEFD TO THE HI-SPEED PRINTER 
CGI IS DDEFD TO 10 UNIT S 
EG1 IS DDEFD TO 10 UNIT 9 
PFRUZ1 IS DDEFD TO 10 UNIT 10 
PFRUY1 IS DDEFD TO 10 UNIT 11 
PFRWZ1 IS DDEFD TO 10 UNIT 12 
PFRWY1 IS DDEFD TO 10 UNIT 13 
CFR1 IS DDEFD TO 10 UNIT 1* 

PP1 IS DDEFD TO 10 UNIT 15 
SSI IS DDEFD TO 10 UNIT 16 
FFG1 IS DDEFD TO 10 UNIT 17 
RUN 1 IS LIBRARY GRAPHICS 
RUN2 IS LIBRARY AESLIB 


AESOP 

DO YOU WISH TO MAKE PLOTS , Y OR N? 

Y 

ENTER THE PLOT NAME - 8 ALPHANUMERIC CHARACTERS 
LUCYTEST 

DO YOU WISH TO MAKE ONLINE PLOTS, Y OR N? 

M 

ENTER TODAY'S DATE (LESS THAN OR EQUAL TO 20 CHARACTERS) 
DEC 2, 1982 

ENTER THE PARAMETER YOU USED FOR AESRUN 
1 

EXTENDED TERMINAL OUTPUT? 


User calls the AESOP program 


READ IN N1 FROM STORAGE? 

DDEF 21 TO THE N1 DATASET *** 

PROCEEDING: PAUSE . PROGRAM WILL PROCEED IF RUN COMMAND ISSUED. 


User requests maximum amount of data to be displayed at 
User elects to have the IFN function number string read 


terminal 
in from storage 


DDEF FT21F0GI»VS»DEC80 User datadefs "dec80", the dataset containing the NAMELIST N1 for the 

test case, to unit 21. "GO” causes program to continue 
GO 

PROCEEDING(PCS) : EXECUTION CONTINUES AT CHCRWC . (X ' 4CD4 * ) 

INI 
IFN = 


101 

201 

401 

402 

403 

801 

802 

301 

803 

804 

805 

806 

807 

808 

809 

810 

302 

811 

303 

812 

813 

814 

815 

816 

817 

818 

819 

820 

501 

502 

503 

504 

505 

506 

507 

508 

509 

510 

511 

512 

513 

514 

515 

516 

517 

518 

519 

520 

521 

522 

523 

524 

525 

601 

602 

603 

604 

701 

702 

703 

704 

705 

210 

404 

801 

301 

809 

819 

405 

999 

0 






liEND 


FUNCTION 101 
CHANGES TO REFS 

DISPLAY REFS BEFORE MAKING CHANGES? 
Y 

4 REFS 

TSFTR= 1.0 
DT= 0.50D-01 
FI = 0.10D-01 
DELF= 0.20D-01 
ZERMAX= 100.0 
AMPSP = 5X1.0 
AMPSR- 5X1.0 
AMPICX= 50X1.0 
I F= 49 
ISP ACE= 1 
IOUT = 1 
IMEAS= 1 
JINC= 1 
JIND= 1 
ITRMX= 100 
NCURV= 2 
L INLOG= 3 
MSPY= 250X1 
MSPYSP= 25X1 
MSPU= 25X1 
M$ROLY= 250X1 


Function 101 requested to demonstrate procedure for changing 
a "reference” parameter (in NAMELIST "REFS’) 
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MSROLX= 250*1 
MICCL Y= 2500*1 
MICCLX= 2500*1 
MICCLU= 250*1 
MICOLY= 2500*1 
MICOLX= 2500*1 
1 END 

ENTER CHANGES TO NAMELIST REFS (TSFTR, DT, FI, DELF, ZERMAX, AMPSP, AMPSR, 

AMPICX, IF, ISPACE, IOUT, IHEAS, JINC, JIND, ITRNX, NCURV, LINLOG, MSPY, 

MSPYSP , MSPU, MSROLY , MSROLX, MICCLY, MICCLX, MICCLU, MICOLY, MICOLX) 

iREFS FS ZERMAX=500 - 0d0 * END Parameter ZERMAX increased from 100 to 500 

TSFTR= 1.0 
DT= 0.50D-01 
FI= 0.10D-01 
DEL F= 0. 200-01 
ZERMAX= 500.0 
AMPSP = 5*1.0 
AMPSR= 5*1.0 
AMPICX= 50*1.0 
IF= 49 
ISPACE= 1 
I OUT = 1 
IMEAS= 1 
JINC= 1 
JIND= 1 
ITRMX= 100 
NCUR V= 2 
L INI OG= 3 
MSPY= 250*1 
MSPYSP= 25*1 
MSPU= 25*1 
MSR0L Y= 250*1 
M$R0LX= 250*1 
MICCL Y = 2500*1 
MICCLX= 2500*1 
MICCLU= 250*1 
MIC0LY= 2500*1 
MIC0LX= 2500*1 
SEND 

ARE THERE ANYMORE CHANGES, Y OR N? 

N 


FUNCTION 201 

PUT IN THE 3RD ORDER TEST CASE MATRICES BY USING SUBROUTINE MATIN 

****** INPUT MATRICES ****** Function that forms all required matrices for running the test case 

A = 


12 3 


1 -o.ioood oo i.ooo 

2 0.0000 0.0000 

3 0.0000 -1.000 


0.0000 
1 .000 

-0.2000D-Q2 


1 2 


1 0.0000 0.0000 

2 0.0000 1.000 

3 1.000 0.0000 


D= 


1 2 


1 0.0000 0.0000 

2 0.0000 1.000 

3 1.000 0.0000 


12 3 


1 1.000 0.0000 

2 0.0000 0.0000 


0.0000 

1.000 


H = 


1 2 3 


1 1.000 0.0000 


0.0000 


66 



DOUT = 


1 


1 0.0000 

2 0.0000 


CSP = 


1 


1 1.000 
2 0.0000 


QC = 


1 


1 20.00 
2 0.0000 
3 0.0000 


NN = 


1 


1 1.000 

2 0.0000 
3 2.000 


PCINV- 


1 


1 1.000 
2 0.0000 


QQ = 


1 


1 0.0000 
2 0.0000 
3 0.0000 


2 


0.0000 

1.000 


2 


0.0000 

0.0000 


2 


0.0000 
1 .000 
0.0000 


2 


0.0000 
3. 000 
0.0000 


2 


0.0000 

0.1000D 00 


2 


0.0000 
2 . 000 
0.0000 


RRINV= 


1 


1 1.000 


t REFS 

TSFTR= 1.0 
DT= 0.50D-01 
FI = 0.10D-01 
DEL F = 0.20D-01 
ZERMAX= 500.0 
AMPSP= 5^1.0 
AMPSR- 5X1.0 
AMPICX= 50X1.0 
IF= 4 9 
ISPACE" 1 
IOUT= 1 
IMEAS= 1 
JINC= 1 
JIND= 1 
I T RMX= 100 
NCURV= 2 
LINLOG- 3 
MSPY= 250X1 
MSPYSP= 25X1 
MSPU= 25X1 
MSROLY= 250X1 
MSROLX= 250X1 
MICCLY= 2500X1 
MICCLX= 2500X1 
MICCLU= 250X1 
MICOLY= 2500X1 
MICOLX= 2500X1 
tEND 


3 


0.0000 

2.000 


3 


0.0000 

0.0000 

10.00 


3 


0 . 0000 
0.0000 
20 .00 
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FUNCTION 401 

OPEN LOOP SYSTEM EIGENVALUES 

XXX* MATRIX IS FOUND TO BE REDUCIBLE XXXX 


The following three functions perform analyses on the 
open-loop system: 


EIGENVALUES 


NAT FREQ (HZ) 

0.1592 
0 . 1592D-01 


ZETA 

0 . 1000D-02 
1.000 


FUNCTION 402 

OPEN LOOP EIGENVECTORS AND MODE SHAPES 
MODIFIED EIGENVECTOR MATRIX OF A 



1 

2 

3 


1 

-1.087 

-0.8933 

1.000 


2 

-1 .001 

0 . 9990 

0.00 00 


3 

1 .000 

1.000 

0 . 00 00 


THE 

MATRIX OF 

MODE SHAPES 

IN MAG. AND ANGL E( DEG . ) FORM 



* 

For 

complex pair at 


' 1 

2 1 

3 


I 

0 . 9951 

-174.4 

1.000 


2 

1.0000 

-90.06 

0.0000 


3 

1.000 

A 

0 . 000 0 

i 

0.0000 



L Magnitudes ■ 

Phase angles. 

deg 


FUNCTION 403 

OPEN LOOP SYSTEM ANALYSIS CALCULATIONS 
SYSTEM CONTROL ABILITY 


1 


2 


1 0.5000 0.500olpor complex eigenvalue pair 

2 0.5730D-01 90.00/ 

3 0.9903 -0.9705D-01 


SYSTEM OBSERVABILITY FOR 
1 2 

1 0.9951 174.4 

SYSTEM OBSERVABILITY FOR 
1 2 

1 0.9951 174.4 

2 1.000 0.0000 


(A AND H) 

3 

1.000 

(A AND C) 

3 

1.000 

0.0000 


RESIDUES FOR (A, B, H) SYSTEM 
RESIDUES FOR (A, B, C) SYSTEM 


FUNCTION SOI 

DESIGN A LINEAR QUADRATIC REGULATOR 


Beginning of Riccati equation solution for LQR problem 


EIGENVALUES OF REGULATOR OR FILTER IN FN/ZETA FORM 
NAT FREQ (HZ) ZETA 

0.2382 0.7116 

0.4519 1.000 

55 = 


12 3 


1 19.24 

2 10.42 

3 1.301 


10.42 1.301 
9.203 1.819 
1.819 1.647 
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KC 


12 3 


1 2.301 1.819 3.647 

2 1.042 1.220 0.1819 


FUNCTION 802 

STORE OPTIMAL CONTROL GAINS (KC) ON UNIT 08. 


FUNCTION 301 
FORM A-BKC MATRIX 


FUNCTION 803 

EIGENVALUES OF SYSTEM WITH STATE FEEDBACK 


EIGENVALUES 


NAT FREQ (HZ) ZETA 

0.2382 0.7116 

0.4519 1.000 


FUNCTION 804 

EIGENVECTORS AND MODE SHAPES WITH STATE FEEDBACK 
MODIFIED EIGENVECTOR MATRIX OF AMBKC 


12 3 


1 -0.9899 0.4246D-01 0.1494 

2 1.000 1.000 -0.4092 

3 0.2144 -1.042 1.000 


THE MATRIX OF MODE SHAPES IN MAG. AND ANGLE(DEG. ) FORM 


12 3 


1 0.7006 -132.5 0.1494 

2 1.000 0.0000 -0.4092 

3 0.7519 123.4 1.000 

* - Beginning of error checks for the Riccati solution matrix 


FUNCTION 805 

POSITIVE DEFINITENESS CHECK OF CONTROL RICCATI SOLUTION MATRIX, SS 


EIGENVALUES 


NAT FREQ (HZ) ZETA 

0 ! 4929 -1 ' 0 0 0> — Negative indicates that eigenvalues all have positive real parts; 

4.133 -1.000) therefore SS is positive-definite 


FUNCTION 806 

SYMMETRY CHECK OF CONTROL RICCATI SOLUTION MATRIX, SS 

MAX. SYMMETRY ERROR IN SS = -3.2429E-15 

AVG. ABSOLUTE SYMMETRY ERROR IN SS = 1.71S4E-15 

SYMMETRY ERROR IN SS = 


1 2 3 


1 0.0000 0.0000 0.0000 

2 -0. 4475 D- 15 0.0000 0.0000 

3 -0 . 3243D-14 -0.1465D-14 0.0000 


FUNCTION 807 

RESIDUAL ERROR CHECK OF SS 

MAX. RESIDUAL, R( 1, 2>= -0.2842D-13 
TRACE OF RESIDUAL = 

-0.2977D-13 
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RESIDUAL ERROR MATRIX FOR CONTROL RICCATI EQUATION 


12 3 


1 -0.1421D-13 -0.2842D-13 -0.1577D-13 

2 -0 . 1754D-13 -0.1622D-13 -0.8438D-14 

3 -0 . 8882D-15 0.4441D-15 0.6661D-15 


The residual remaining upon substitution of solution 
matrix SS back into the original Riccati equation 


FUNCTION 808 

STORE CONTROL RICCATI SOLUTION MATRIX (SS) ON UNIT 16 

Beginning of Riccati equation solution for Kalman filter problem 


FUNCTION 809 

DESIGN A KALMAN FILTER 


EIGENVALUES OF REGULATOR OR FILTER IN FN/2ETA FORM 
NAT FREQ (HZ) ZETA 


0.2202 1.000 

0.2862 0.4<t51 


1 2 3 


1 2.882 4.442 

2 4.442 11.76 

3 1.482 8.865 


1.482 

8.865 

18.37 


KE = 


1 


1 2.882 

2 4.442 

3 1.482 


FUNCTION 810 

STORE KALMAN FILTER GAINS (KE) ON UNIT 09 


FUNCTION 302 

FORM A-BKC-KEH MATRIX 


FUNCTION 811 

EIGENVALUES OF OPTIMAL CONTROLLER A-BKC-KEH 


EIGENVALUES 


NAT FREQ (HZ) ZETA 


0.6057 1.000 

0.5366 0.6000 


FUNCTION 303 

FORM ATOT , CTOT, DTOT , KCTOT, AND HTOT MATRICES 

FOR OPTIMAL CONTROL SYSTEM WITH KALMAN FILTER IN FEEDBACK LOOP 

ATOT = 



1 

2 

3 

4 

5 

6 

1 

-0.1000D 00 

1.000 

0.0000 

0.0000 

0.0000 

0.0000 

2 

0.0000 

0= 0000 

1.000 

-1.042 

-1.220 

-0.1819 

3 

0.0000 

” 1.000 

-0.2000D-02 

-2.301 

-1.819 

-3.647 

4 

2.882 

0.0000 

0.0000 

-2.982 

1.000 

0.0000 

5 

4.442 

0.0000 

0.0000 

-5.484 

-1.220 

0.8181 

6 

1.482 

0.0000 

0.0000 

-3.783 

-2.819 

-3.649 

CTOT 

= 







X 

2 

3 

4 

5 

6 

1 

1.000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

2 

0.0000 

0.0000 

1.000 

-1.042 

-1.220 

-0.1819 
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DTOT 


1 2 


1 0.0000 0.0000 

2 0.0000 1.000 

3 1.000 0.0000 

4 0.0000 0.0000 

5 0.0000 0.0000 

6 0.0000 0.0000 


KCTOT □ 



1 

2 

3 4 

5 

6 

1 

0.0000 

0.0000 

0.0000 -2.301 

-1.819 

-3.647 

2 

0.0000 

0.0000 

0.0000 -1.042 

-1.220 

-0.1819 

HTOT 

= 






1 

2 

3 4 

5 

6 

1 

1.000 

0.0000 

0.0000 0.0000 

0.0000 

0.0000 

FUNCTION 812 
EIGENVALUES OF 

CONTROL 

SYSTEM WITH A FILTER IN THE 

FEEDBACK 

LOOP 


EIGENVALUES 


NAT FREQ (HZ) ZETA 


0.2202 

0.2382 

0.2862 

0.4519 


1 . 000 
0.7116 
0.4451 
1.000 




These are 
These are 


the eigenvalues of A - KE 
the eigenvalues of A - B • 


* H 
KC 


FUNCTION 813 ’♦—Beginning of error checks on Kalman filter error covariance matrix 

POSITIVE DEFINITENESS CHECK OF ERROR COVARIANCE MATRIX, PP 


EIGENVALUES 


NAT FREQ (HZ) ZETA 

0.1103 -1.000) 

1.138 -1.000} Matrix PP is positive-definite 

4.006 -1.000) 


FUNCTION 814 

SYMMETRY CHECK OF ERROR COVARIANCE MATRIX, PP 

MAX. SYMMETRY ERROR IN PP = -2.9962E-16 

AVG. ABSOLUTE SYMMETRY ERROR IN PP = 1.6667E-16 

SYMMETRY ERROR IN PP = 


12 3 


1 0.0000 0.0000 0.0000 
2 0.0000 0.0000 0.0000 
3 -0.2996D-15 0.2004D-15 0.0000 


FUNCTION 815 

RESIDUAL ERROR CHECK OF PP 

MAX. RESIDUAL, R( 2, 2)= -0.1155D-13 
TRACE OF RESIDUAL = 

-0.2265D-13 

RESIDUAL ERROR MATRIX FOR ESTIMATION RICCATI EQUATION 


12 3 


1 -0 . 4441D-15 -0.3775D-14 -0.5958D-14 

2 -0.4219D-14 -0 . 1155D-13 -0.1066D-13 

3 -0.2849D-14 -0.1Q87D-13 -0.1066D-13 



FUNCTION 816 

STORE ERROR COVARIANCE MATRIX (PP) ON UNIT 15 


FUNCTION 817 

COMPUTE COVARIANCE MATRICES FOR LINEAR QUADRATIC REGULATOR 
WITH KALMAN FILTER IN FEEDBACK LOOP 


XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

COVARIANCE MATRICES 


xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


UU, CONTROL COVARIANCE MATRIX 


1 2 


1 41.42 20.81 

2 20.81 25.13 


XX, STATE COVARIANCE MATRIX 


1 2 3 


1 20.99 2.099 

2 2.099 21.56 

3 -6.849 8.428 


-6.849 

8.428 

24.26 


YY, OUTPUT COVARIANCE MATRIX 


1 2 


1 20.99 -21.35 

2 -21.35 65.67 


ZZ, MEASUREMENT COVARIANCE MATRIX 


1 


1 20.99 


FUNCTION 818 

LYAPUNOV ERROR CHECK FOR FUNCTION 817 * Computes 

THE TRACE OF THE RESIDUAL 2 0.65148D-12 
NORMALIZED DIAGONAL ELEMENTS OF THE ERROR MATRIX = 

-0 . 6497 0D-15 0 . 6 3 0 51D- 15-0 . 4 0 6 17D-14 
THE TRACE OF THE ERR0R=-0 . 98575D-13 
THE TRACE OF THE COVARIANCE 2 66.806 
TR(ERROR)/TR(COV. ) =-0.147550-14 


error incurred in calculating state covariance matrix XX 


FUNCTION 819 

FORM FEED FORWARD MATRIX FOR NON-ZERO SET POINT CONTROL 
KFF 2 


1 2 


1 2.583 1.824 

2 1.164 -0.4090 


FUNCTION 820 

STORE FEED FORWARD MATRIX (KFF) ON UNIT 17 


Beginning of frequency response and Bode plot computation 


FUNCTION 501 

OPEN LOOP FREQUENCY RESPONSE OF MEASUREMENT 1 TO CONTROL INPUT 1 


NUMERATOR 

COEFFICIENTS 

DENOMINATOR COEFFICIENTS 

0.00000 

xsxx 3 

1.0000 

xsxx 

3 

-0.12089D- 

-16XSXX 2 

0.10200 

xsxx 

2 

0.00000 

xsxx 1 

1.0002 

xsxx 

1 

1.0000 

xsxx 0 

0.10000D 

OOXSxx 

0 
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II III 


I I II I I 


III ■ 1 1 


FUNCTION 502 

PLOT OPEN LOOP FREQ. RESPONSE OF MEASUREMENT TO CONTROL INPUT 


GRAPHICS DEVICE NOT DEFINED BY DDEF . 

ENTER UNIT NAME. DEFAULT TO CANCEL. 

2ETA12 

PLOT IDENTIFICATION = PELUCYXX 02 DEC 1982 9 08:06U8 


User specifies ZETA plotter to be used for plotting 


FUNCTION 503 

STORE OPEN LOOP FREQ. RESPONSE OF MEASUREMENT TO CONTROL INPUT ON UNIT 10 


FUNCTION 504 

OPEN LOOP FREQUENCY RESPONSE OF OUTPUT 1 TO CONTROL INPUT I 


NUMERATOR 

COEFFICIENTS 

DENOMINATOR COEFFICIENTS 

0.00000 

XSXX 

3 

1.0000 xsxx 

3 

-0.12089D 

- 1 6 x s x x 

2 

0.10200 xsxx 

2 

0.00000 

xsxx 

1 

1.0002 xsxx 

1 

1.0000 

xsxx 

0 

0.10000D 00XSXX 

0 

FUNCTION 

505 





PLOT OPEN LOOP FREQ. RESPONSE OF OUTPUT TO CONTROL INPUT 


FUNCTION 506 

STORE OPEN LOOP FREQ. RESP0N5E OF OUTPUT TO CONTROL INPUT ON UNIT 11 


FUNCTION 507 

OPEN LOOP FREQUENCY RESPONSE OF MEASUREMENT 1 TO DISTURBANCE INPUT 1 


NUMERATOR 

COEFFICIENTS 

DENOMINATOR COEFFICIENTS 

0 .00000 

xsxx 

3 

1 .0000 

XSXX 

3 

-0. 120S9D- 

16 xsxx 

2 

0.10200 

XSXX 

2 

0.00000 

xsxx 

1 

1 .0002 

XSXX 

1 

1.0000 

xsxx 

0 

0 .10000D 

OOXSxx 

0 

FUNCTION 

508 






PLOT OPEN LOOP FREQ. RESPONSE OF MEASUREMENT TO DISTURBANCE INPUT 


FUNCTION 509 

STORE OPEN LOOP FREQ. RESPONSE OF MEASUREMENT TO DISTURBANCE INPUT ON UNIT 12 


FUNCTION 510 

OPEN LOOP FREQUENCY RESPONSE OF OUTPUT 1 TO DISTURBANCE INPUT 1 


NUMERATOR 

COEFFICIENTS 

DENOMINATOR COEFFICIENTS 

0.00000 

xsxx 

3 

1.000 0 

xsxx 

3 

-0 .12089D- 

•16XSXX 

2 

0.10200 

xsxx 

2 

0.00000 

xsxx 

1 

1 .0002 

xsxx 

1 

1 .000 0 

xsxx 

0 

0. 10000D 

OOxsxx 

0 

FUNCTION 

511 






PLOT OPEN LOOP FREQ. RESPONSE OF OUTPUT TO DISTURBANCE INPUT 


FUNCTION 512 

STORE OPEN LOOP FREQ. RESPONSE OF OUTPUT TO DISTURBANCE INPUT ON UNIT 13 


FUNCTION 513 

CLOSED LOOP FREQUENCY RESPONSE OF OUTPUT 1 TO DISTURBANCE INPUT 1 
FOR STATE FEEDBACK 


NUMERATOR COEFFICIENTS 

0.00000 XSxx 3 

-0. 22204 D-15XSXX 2 
-0 .22204D-15XSX* 1 

0.81810 xsxx 0 


DENOMINATOR COEFFICIENTS 
1.0000 #Sxx 3 

4.9693 *SK* 2 

8.2882 xsxx l 

6.3604 xsxx 0 


FUNCTION 514 

PLOT CLOSED LOOP FREQ. RESPONSE OF OUTPUT TO DISTURBANCE INPUT 
FOR STATE FEEDBACK 

(PLUS THE CORRESPONDING OPEN LOOP RESPONSE, IF DESIRED) 


FUNCTION 515 

CLOSED LOOP FREQUENCY RESPONSE OF CONTROL 1 TO DISTURBANCE INPUT I 
FOR STATE FEEDBACK 


NUMERATOR COEFFICIENTS 


0.00000 
-3.6470 
-4. 3034 
-6.2763 


xsxx 

xsxx 

xsxx 

xsxx 


DENOMINATOR COEFFICIENTS 
1.0000 xsxx 3 

4.9693 xsxx 2 

8.2882 xsxx l 

6.3604 xsxx 0 
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FUNCTION 516 

PLOT CLOSED LOOP FREQ. RESPONSE OF CONTROl TO DISTURBANCE INPUT 
FOR STATE FEEDBACK 


FUNCTION 517 

CLOSED LOOP FREQUENCY RESPONSE OF MEASUREMENT 1 TO DISTURBANCE INPUT 1 
FOR CONTROL SYSTEM WITH A FILTER IN THE FEEDBACK LOOP 

NUMERATOR COEFFICIENTS DENOMINATOR COEFFICIENTS 


0.00000 

XSXX 

6 

1 . 0000 

xsxx 

6 

0.00000 

xsxx 

5 

7.9536 

xsxx 

5 

0.00000 

XSXX 

4 

28.566 

xsxx 

4 

1.0000 

xsxx 

3 

62.639 

xsxx 

3 

7.8516 

xsxx 

2 

86.362 

xsxx 

2 

26.764 

XSXX 

1 

71.725 

xsxx 

1 

43.263 

xsxx 

0 

28.452 

xsxx 

0 


FUNCTION 518 

PLOT CLOSED LOOP FREQ. RESPONSE OF MEASUREMENT TO DISTURBANCE INPUT 

FOR CONTROL SYSTEM WITH FILTER IN FEEDBACK LOOP 

(PLUS THE CORRESPONDING OPEN LOOP RESPONSE, IF DESIRED) 


FUNCTION 519 

CLOSED LOOP FREQUENCY RESPONSE OF OUTPUT 1 TO DISTURBANCE INPUT 1 
FOR CONTROL SYSTEM WITH A FILTER IN THE FEEDBACK LOOP 


NUMERATOR 

COEFFICIENTS 

DENOMINATOR COEFFICIENTS 

0.00000 

xsxx 

6 

1 .0000 

xsxx 

6 

0.00000 

xsxx 

5 

7.9536 

xsxx 

5 

0.00000 

xsxx 

4 

28.566 

xsxx 

4 

1.0000 

xsxx 

3 

62.639 

xsxx 

3 

7.8516 

xsxx 

2 

86 .362 

xsxx 

2 

26.764 

xsxx 

1 

71.725 

xsxx 

1 

43.263 

xsxx 

0 

28.452 

xsxx 

0 


FUNCTION 520 

PLOT CLOSED LOOP FREQ. RESPONSE OF OUTPUT TO DISTURBANCE INPUT 

FOR CONTROL SYSTEM WITH FILTER IN FEEDBACK LOOP 

(PLUS THE CORRESPONDING OPEN LOOP RESPONSE, IF DESIRED) 


FUNCTION 521 

CLOSED LOOP FREQUENCY RESPONSE OF CONTROL 1 TO DISTURBANCE INPUT 1 
FOR CONTROL SYSTEM WITH A FILTER IN THE FEEDBACK LOOP 

NUMERATOR COEFFICIENTS DENOMINATOR COEFFICIENTS 


0.00000 

xsxx 

6 

1 . 0000 

xsxx 

6 

-0.44409D- 

-15X5XX 

5 

7.9536 

xsxx 

5 

-0. 13323D- 

-14xsxx 

4 

28.566 

xsxx 

4 

-0 .39968D- 

-14XSXX 

3 

62.639 

xsxx 

3 

-20.117 

xsxx 

2 

86.362 

xsxx 

2 

-6.8315 

xsxx 

1 

71.725 

xsxx 

1 

-24.088 

xsxx 

0 

28.452 

xsxx 

0 

FUNCTION 

522 






PLOT CLOSED LOOP FREQ. RESPONSE OF CONTROL TO DISTURBANCE INPUT 
FOR CONTROL SYSTEM WITH FILTER IN FEEDBACK LOOP 


FUNCTION 523 

FREQUENCY RESPONSE OF CONTROL 1 TO MEASUREMENT 1 


FOR OPTIMAL CONTROLLER 


NUMERATOR 
0.00000 
-20.117 
-6.8315 
-24 . 088 


COEFFICIENTS 
XSXX 3 
XSXX 2 
KS*x 1 
xs xx 0 


DENOMINATOR COEFFICIENTS 
1.0000 XSXX 3 

7.8516 xsxx 2 

26.764 xsxx 1 

43.263 XSXX Q 


FUNCTION 524 

PLOT FREQ. RESPONSE OF CONTROL TO MEASUREMENT 
FOR OPTIMAL CONTROLLER 


FUNCTION 525 

STORE FREQ. RESPONSE OF CONTROL TO MEASUREMENT 
FOR OPTIMAL CONTROLLER ON UNIT 14 


Beginning of transient response calculations and plots 


FUNCTION 601 

OBTAIN AND PLOT SELECTED OPEN LOOP STEP RESPONSES 

STATE TRANSITION MATRIX FOR OPEN LOOP SYSTEM FOR TIME STEP =Q.50Q0D-0I 
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1 


2 


3 


1 0.9950 

2 0.0000 

3 0.0000 


0.4985D-01 
0.9938 
-0 . 4998D-01 


0 . 1248D-02 
0.4993D-01 
0 . 9987 


FORCED RESPONSE MATRIX OF OPEN LOOP SYSTEM FOR TIME STEP =0.500QD-01 


1 2 3 


1 

2 

3 


0 . 4988D-01 
0.0000 
0.0000 


0 . 1248D-02 
0.499SD-01 
-0.I250D-02 


0.2080D-04 
0.1250D-02 
0 . 4 9 98D- 0 1 


FUNCTION 602 

OBTAIN AND PLOT SELECTED INITIAL CONDITION RESPONSES FOR THE OPEN LOOP SYSTEM 
STATE TRANSITION MATRIX FOR OPEN LOOP SYSTEM FOR TIME STEP =0.50000-01 


12 3 


1 0.9950 

2 0.0000 

3 0.0000 


0.4985D-01 0 . 1248D-Q2 

0.9988 Q.4998D-Q1 

-0.4998D-01 0.9987 


FORCED RESPONSE MATRIX OF OPEN LOOP SYSTEM FOR TIME STEP =0.5000D-01 


12 3 


1 

2 

3 


0 .4988D-01 
0.0000 
0.0000 


Q.1248D-02 
0 .4998D-01 
-0 . 125 0 D- 02 


0 .20S0D-04 
0 .1250D-02 
0 .499SD-01 


FUNCTION 603 

OBTAIN AND PLOT SELECTED INITIAL CONDITION RESPONSES 
FOR THE CLOSED LOOP LINEAR REGULATOR 

STATE TRANSITION MATRIX OF LINEAR REGULATOR FOR TIME STEP 0.5000D-01 


12 3 


1 0.9937 0 . 4832D-Q1 

2 -0.5251D-01 0.9369 

3 -0.1014 -0.1273 


0. 9413D-03 
0.3619D-01 
0.8307 


FORCED RESPONSE MATRIX OF LINEAR REGULATOR FOR TIME STEP 0.5000D-01 


12 3 


1 0 . 4985D-01 0 . 1222D-Q2 

2 -0.1310D-02 0 . 4844D- 0 1 

3 -0.2645D-02 -Q.3294D-02 


0 .1602D-04 
0. 9429D-03 
0 . 4566D-01 


FUNCTION 604 

OBTAIN AND PLOT SELECTED STEP RESPONSES 
FOR THE NON-ZERO SET POINT LINEAR REGULATOR 

STATE TRANSITION MATRIX OF LINEAR REGULATOR FOR TIME STEP 0.5000D-01 


12 3 


1 0.9937 0.4S32D-01 

2 -0 . 5251D-01 0.9369 

3 -0.1014 -0.1273 


0 . 9413D-03 
0 . 3619D-0 1 
0.8307 


FORCED RESPONSE MATRIX OF LINEAR REGULATOR FOR TIME STEP 0.500QD-01 


12 3 


1 0 . 4985D-01 0 . 1222D-02 

2 -0.1310D-02 0.4844D-01 

3 -0.2645D-02 -0.3294D-02 


0 . 1602D-04 
0 . 9429D-03 
0 . 4566D-01 
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-«■ Beginning of transfer function gain and numerator zeros calculation 

FUNCTION 701 

GAIN AND ZEROES OF OPEN LOOP TRANSFER FUNCTION 
RELATING MEASUREMENT 1 TO CONTROL INPUT 1 

GAIN = 1.00000 

NUMBER OF ZEROES = 0 

REAL PARTS OF NUMERATOR ZEROES 


1 


1 0.0000 


IMAGINARY PARTS OF NUMERATOR ZEROES 


1 


1 0.0000 


FUNCTION 702 > 

GAIN AND ZEROES OF OPEN LOOP TRANSFER FUNCTION 
RELATING OUTPUT 1 TO CONTROL INPUT 1 

GAIN = 1.00000 


NUMBER OF ZEROES = 0 

REAL PARTS OF NUMERATOR ZEROES 

1 

1 0.0000 

IMAGINARY PARTS OF NUMERATOR ZEROES 
1 

1 0.0000 


FUNCTION 703 

GAIN AND ZEROES OF OPEN LOOP TRANSFER FUNCTION 
RELATING MEASUREMENT 1 TO DISTURBANCE INPUT 1 

GAIN = 1.00000 

NUMBER OF ZEROES = 0 

REAL PARTS OF NUMERATOR ZEROES 


1 


1 0.0000 


IMAGINARY PARTS OF NUMERATOR ZEROES 


1 


1 0.0000 


FUNCTION 704 

GAIN AND ZEROES OF OPEN LOOP TRANSFER FUNCTION 
RELATING OUTPUT 1 TO DISTURBANCE INPUT 1 

GAIN = 1.00000 

NUMBER OF ZEROES = 0 

REAL PARTS OF NUMERATOR ZEROES 


1 

1 0.0000 

IMAGINARY PARTS OF NUMERATOR ZEROES 


1 


1 0.0000 
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FUNCTION 705 

GAIN AND ZEROES OF OPTIMAL CONTROLLER TRANSFER FUNCTION 
RELATING CONTROL 1 TO MEASUREMENT 1 

GAIN = 20.1170 

NUMBER OF ZEROES = 2 

REAL PARTS OF NUMERATOR ZEROES 


1 


1 -0.1699 

2 -0.1699 

IMAGINARY PARTS OF NUMERATOR ZEROES 


1 


1 1.081 

2 -1.081 Beginning of series of functions that demonstrate (1) normalizing of input matrices, 

(2) normalizing of control weighting matrices, (3) recomputing of KC, KE, and KFF 
matrices, and (4) unnormalizing KC, KE, KFF, and PP to show that results are 
FUNCTION 210 identical to prior calculations 

CHANGE MATRICES A, B, C, D, H, DOUT, CSP, QC, NN, PCINV, QQ, OR RRINV AND ANY 
OR ALL DIMENSIONS BY READING CHANGES IN AT TERMINAL USING NAMELIST * MATDAT ' 

DISPLAY MATDAT BEFORE MAKING CHANGES? 

Y 

AMATDAT 

A= -0.10D0, 49*0.0, 1.0, 0.0, -1.0, 48*0.0, 1.0, -0.20D-02, 2397*0.0 

B= 2*0.0, 1.0, 48*0.0, 1.0, 198*0.0 

C= 1.0, 100*0.0, 1.0, 2398*0.0 

D= 2*0.0, 1.0, 48*0.0, 1.0, 698*0.0 

H= 1.0, 249*0.0 

DOUT = 51*0.0, 1.0, 198*0.0 

CSP= 1.0, 10*0.0, 2.0, 238*0.0 

QC = 20.0, 50*0.0, 1.0, 50*0.0, 10.0, 2397*0.0 

NN= 1.0, 0.0, 2.0, 48*0.0, 3.0, 198*0.0 

PC I N V = 1.0, 5*0.0, 0.10D0, 18*0.0 

QQ = 51*0.0, 2.0, 50*0.0, 20.0, 2397*0.0 

RR I N V= 1.0, 24*0.0 

N= 3 

NM= 1 

NC= 2 

ND= 2 

NO= 2 


tEND 

ENTER MATDAT CHANGES AS tMATDAT A<1,1)= 
X MATDAT 

QC C 1 , 1 ) = 500. 0D0 
QC ( 2 , 2 ) = 9.0D0 
QC ( 3 , 3 ) = 4 0 . 0 DO 
PCINV ( 1 , 1 ) = 1.5625D-2 
PCINV(2,2) = 1.2345679D-3 
NN ( 1 , 1 ) = 4 0 . 0 DO 
NN ( 3 , 1 ) = 32 . 0 DO 
N N ( 2 , 2 ) = 8 1 . 0 DO 
t END 


, ETC. tEND 

• New values for weighting elements calculated off-line so 

that they correspond to normalized versions of those 
stored above as MATDAT 


iMATDAT 

A- -0.10D0, 49*0.0, 1.0, 0.0, -1.0, 48*0.0, 1.0, -0.20D-02, 

B= 2*0.0, 1.0, 48*0.0, 1.0, 198*0.0 

C= 1.0, 100*0.0, 1.0, 2398*0.0 

D= 2*0.0, 1.0, 48*0.0, 1.0, 698*0.0 

H= 1.0, 249*0.0 

DOUT = 51*0.0, 1.0, 198*0.0 

CSP= 1.0, 10*0.0, 2.0, 238*0.0 

QC= 500.0, 50*0.0, 9.0, 50*0.0, 40.0, 2397*0.0 

NN= 40.0, 0.0, 32.0, 48*0.0, 81.0, 198*0.0 

PCINV= 0 . 156250D-01, 5*0.0, 0 . 1234567 90D-02 , 18*0.0 

QQ= 51*0.0, 2.0, 50*0.0, 20.0, 2397*0.0 

RRINV= 1.0, 24*0.0 

N= 3 

NM- 1 

NC= 2 

ND= 2 

NO= 2 


AEND 

ARE THERE ANYMORE CHANGES, Y OR N? 
N 


2397*0.0 


FUNCTION 404 

NORMALIZE SYSTEM MATRICES; 

NORMALIZING FACTORS ARE READ IN FROM UNIT 34 IF NOT PREVIOUSLY SET 

IF HOT ALREADY DONE, DDEF DATASET CONTAINING NAMELIST * NRMS * TO UNIT 34 
PROCEEDING ' PAUSE . PROGRAM WILL PROCEED IF _RUN COMMAND ISSUED. 

DDEF FT34F001,V5, ANRMS 

GO 

PROCEEDING* PCS >*• EXECUTION CONTINUES AT CHCRWC . (X* 4CD4* ) 


NORMALIZING FACTORS ARE 
SNRMS 

SCX= 5.0, 3.0, 2.0, 47*0.0 
SCU= 8.0, 9.0, 3*0.0 



SCY= 4.0, 7.0, 10.0, 47X0.0 
SCZ= 6.0, 4X0.0 
SCY5P = 12.0, 11.0, 3X0.0 
& END 

NORMALIZED SYSTEM MATRICES 


THE 

A MATRIX 

IS 




1 


2 

3 

1 

-0 .1000D 

00 

0.6000 

0.0000 

2 

0 .0000 


0.0000 

0.6667 

3 

0 . 0000 


-1.500 

-0 .200 0D- 

THE 

B MATRIX 

IS 




1 


2 


1 

0.0000 


0.0000 


2 

0.0000 


3.000 


3 

4.000 


0.0000 


THE 

C MATRIX 

IS 




1 


2 

3 

1 

1.250 


0.0000 

0.0000 

2 

0.0000 


0.0000 

0.2857 

THE 

H MATRIX 

IS 




1 


2 

3 

1 

0.8333 


0.0000 

0.0000 

THE 

QQ MATRIX 

IS 




1 


2 

3 

1 

0.0000 


0.0000 

0.0000 

2 

0.0000 


0.2222 

0.0000 

3 

0.0000 


0.0000 

5.000 


THE RRINV MATRIX IS 
1 

1 36.00 
THE D MATRIX IS 

1 2 

1 0.0000 0.0000 

2 0.0000 0.3333 

3 0.5000 0.0000 

THE DOUT MATRIX IS 

1 2 

1 0.0000 0.0000 

2 0.0000 1.286 

THE CSP MATRIX IS 

12 3 

1 0.4167 0.0000 0.0000 

2 0.0000 0.0000 0.3636 

FUNCTION 801 

DESIGN A LINEAR QUADRATIC REGULATOR 


78 



EIGENVALUES OF REGULATOR OR FILTER IN FN/ZETA FORM 
NAT FREQ (HZ) ZETA 


SS = 


0.2382 0.71161 

0.4519 1.000/ 


Note that LQR eigenvalues are not changed by process of 
normalization (QC, NN, and PCINV were also normalized) 


12 3 


1 481 . 1 156.3 

2 156.3 82.83 

3 13.01 10.91 


13.01 

10.91 

6.588 


KC = 


12 3 


1 1.438 0.6821 

2 0.5789 0.4068 


0.9117 
0 .40420-01 


FUNCTION 301 
FORM A-BKC MATRIX 


FUNCTION 809 

DESIGN A KALMAN FILTER 


EIGENVALUES OF REGULATOR OR FILTER IN FN/ZETA FORM 


NAT FREQ (HZ) ZETA 

0.2202 1.000 

0.2862 0.4451 


Same Kalman filter eigenvalues as obtained previously 


12 3 


1 0.1153 0.2961 

2 0.2961 1.307 

3 0.1482 1.477 


0.1482 

1.477 

4.591 


KE = 


1 


1 3.459 

2 8.884 

3 4.446 


FUNCTION 819 

FORM FEED FORWARD MATRIX FOR NON-ZERO SET POINT CONTROL 
KFF = 

1 2 

1 3.874 2.509 

2 1.552 -0.4999 


FUNCTION 405 

UNNORMALIZE GAINS AND ERROR COVARIANCE MATRIX; 

NORMALIZING FACTORS ARE READ IN FROM UNIT 34 IF NOT PREVIOUSLY SET 


NORMALIZING FACTORS 
tNRMS 

SCX= 5.0, 3.0, 2.0, 47*0.0 

SCU= 8.0, 9.0, 3*0.0 

SCY- 4.0, 7.0, 10.0, 47*0.0 

SCZ= 6.0, 4*0.0 

SCYSP= 12.0, 11.0, 3*0.0 

IEND 

Check shows that the following matrices are identical to those previously computed: 


UNNORMALIZED SYSTEM MATRICES 
THE KC MATRIX IS 



1 


2 


3 


1 2.301 1.819 3.647 

2 1.042 1.220 0.1819 

THE KE MATRIX IS 
1 

1 2.882 

2 4.442 

3 1.482 

THE KFF MATRIX IS 

1 2 

1 2.583 1.824 

2 1.164 - 0.4090 

THE PP MATRIX IS 


2.882 

4.442 

1.482 


4.442 

11.76 

8.865 


1.482 

8.865 

18.37 


STORE THE HI FOR THIS RUN? if new function number string has been formed, it could have been stored 

Terminated: stop return at this P° int for future use 

PRINT 0UT1 , PRTSP=EDIT .*-User requests output dataset (OUTl in this case) to be printed out 


Test Case II - Interactive Design of a Nonzero-Set -Point Regulator 

Test case II demonstrates how AESOP can be used in an interactive manner to design a feedback 
control system. The same third-order state-variable system model used in test case I is considered, 
except that the outputs of primary interest are now the set -point outputs y sp . They are selected to be 
y sp i = x\ and y sp 2 = x 3 by proper definition of matrix CSP. The following assumptions are made: 

(1) All three states are measurable. 

(2) Both disturbance w and measurement noise v are zero. 

The design problem is to compute feedforward (KFF) and feedback (KC) gain matrices such that the 
resulting closed-loop system meets the following design criteria: 

(1) Each of the two set -point outputs follows a step in its corresponding set point with zero steady- 
state error. (This is insured by the nonzero-set -point regulator structure.) 

(2) Set -point output step reponses are well damped (less than 10 percent overshoot) and settle out 
in less than 5 seconds. 

(3) For a unit step on either set point, excursions of the two control variables must be such that 
|uj|<5 and |ii2j^5. 

The resulting closed-loop system is shown in figure 15. Basically, the design procedure followed is to 
vary performance index weights in an interactive fashion, displaying the step responses of y sp j and 
y sp2 for each candidate design and repeating the process until acceptable transient responses are 
observed. Once transient performance is considered acceptable (criteria 1 and 2 are satisfied), the 
control variable responses are displayed to check that criterion 3 is also met. 

The input data that define the problem are stored in dataset TEST.CASE2, which is shown in 
figure 16. Note that the values of some of the reference parameters in NAMELIST REFS have been 
made different from the default values so as to obtain the desired output variable selection and time 
step for the step responses. Also, unlike test case I the initial values of weighting matrices QC and NN 
are allowed to be zero, and PCINV is set equal to a 2x2 identity matrix. 

To begin the design process the AESOP program is initiated as was done in test case I. When the 
program prompts the user with the message READ IN N1 FROM STORAGE?, the user replies “N.” 
The program then prompts the user to enter function numbers from the terminal 
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0488900 

0489000 

0489100 

0489200 

0489300 

0489400 

0489500 

0489600 

0489700 

0489800 

0489900 


&MATDAT 

N=3,NC=2,NM=1,N0=2,ND=2, 

AC1,1)=-0.1,0.,0.,A(1,2)=1.,0.,-1.,A(1,3)=0.,1. ,“.2D~2, 
BC3,1)=1.,B(2,2)=1.,D(3,1)=1.,D(2,2)=1.,C(1,1)=1.» 
C(2,3)=1.,H<1,1)=1. ,D0UT(2.2)=1.,CSPC1,1)=1. ,CSPC2,3)=1. , 
PCINV<1,1)=1.,PCINV(2,2)=1.,QQ<2,2)=2. , QQC 3 , 3) =20 . , RRINVC 1 , 1 > =1 . 
& END 
&REFS 
DT= . 24 , 

MSPY = 25 0K0,MSPYSP=25*0,MSPU=25#O,MSPYSPU, 1 ) =1 , MSPYSP (2 , 2 ) = 1 
& END 


Figure 16. - Data for test case II, contained in dataset TEST.CASE2. 


ENTER NAMELIST DATA AS ‘&N1 IFN = , , , &END 


and the user responds with 


&nl ifn = 202,203,801 &end 


The user then datadefs dataset TEST.CASE2 to unit 33. 


FUNCTION 202 

READ INPUT DATA — MATRICES AND REFERENCE VALUES DEFINED 
IN NAMELISTS ' MATDAT ' AND 'REFS' 

IF NOT ALREADY DONEt DDEF DATASET CONTAINING NAMELISTS 

'MATDAT' AND 'REFS' TO UNIT 33 

CHCRU410 proceeding: PAUSE 

ddef f t33f 001 t vs ? test • csse2 

2o 

CZAPB030 PROCEEDING ( PCS ) t EXECUTION CONTINUES AT CHCRWC < < X ' 4CD4 ' ) 

DISPLAY INPUT MATRICES? 
n 


The next requested function (203) allows the user to enter desired performance index weights. For 
this test case, weights on the two controls will be fixed at unity and the weight on state 1 will be 
varied. A low value (0.001) is selected initially, and the program proceeds to compute an LQR 
solution. 


FUNCTION 203 

DISPLAY CONPAR BEFORE MAKING CHANGES? 
* 

1CONPAR 
GC= 2500*0 * 0 
NN= 250*0,0 


81 







PC I NV = l.Or 5*0*0* 1*0* 18*0 ♦ 0 
SEND 

ENTER CONTROL WTS GC* NN AND/OR PCINV (NAMELIST CONPAR) 
Seonpsr ac(l*l) = 0*001 Send 
SCONPAR 

QC= 0.10D-02* 2499*0*0 
NN= 250*0*0 

P C I N 0 = 1*0» 5*0 ♦ 0 * 1.0* 18*0*0 
SEND 

ARE THERE ANYMORE CHANGES * Y OR N? 
n 


FUNCTION 801 

DESIGN A LINEAR QUADRATIC REGULATOR 


EIGENVALUES OF REGULATOR OR FILTER IN FN/ZETA FORM 
NAT FREQ <HZ> ZETA 


SS 


0* 1667 D -01 1*000 

0*1593 0*2223 D -01 


1 


3 


1 

2 

3 


0*4883D-02 
0.61550-03 
0 * 48060-02 


0. 6155 D -03 
0.2122D-01 
0 ♦ 3902D-03 


0*48060-02 

0*39020-03 

0*26010-01 


KC = 


1 


3 


1 

2 


0*48060-02 

0*61550-03 


0*39020-03 

0 . 21220-01 


0*26010-01 

0*39020-03 


FUNCTION 0 

TO COMPUTE FURTHER* ENTER NEXT FUNCTION N0S*(I3>* UNE PER LINE 
TO TERMINATE ENTER 999<LAST ENTRY MUST BE A RETURN) 

301 

819 

604 


The user then requests that step responses be plotted for a nonzero-set-point regulator. This 
requires, as prerequisites, forming of A - B-KC and calculation of feedforward gain matrix KFF. 


FUNCTION 301 
FORM A-BKC MATRIX 


FUNCTION 819 

FORM FEED FORWARD MATRIX FOR NON-ZERO SET POINT CONTROL 
KFF « 


1 


2 


1 0*1048 0.2801D-01 

2 0. 27370-02 -0*9996 
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FUNCTION 604 

OBTAIN AND PLOT SELECTED STEP RESPONSES 
FOR THE NON-ZERO SET POINT LINEAR REGULATOR 

STATE TRANSITION MATRIX OF LINEAR REGULATOR FOR TIME STfcP 0*2400 


1 


3 


1 0.9763 0.2343 

2 -0 • 2602D-03 0.9664 

3 -0.1 1 08D-02 -0.2365 


0.2831 D -01 

0.2362 

0.9648 


FORCED RESPONSE 


MATRIX OF LINEAR REGULATOR FOR TIME STEP 0.2400 


1 


2 


3 


1 0.2371 0.2839D-01 

2 -0 ♦ 2841D-04 0.2371 

3 -0.1349D-03 -0.2857D-01 


0 . 2276 D -02 
0 » 2854D-0 1 
0.2369 


Next, before plots can be obtained, the user must indicate on which device the plots are to be 
generated by entering the user’s terminal number (LA001), indicating that plots are to be displayed 
not on an off-line device but at the user’s terminal. 


GRINT100 GRAPHICS DEVICE NOT DEFINED BY DDEF . 

ENTER UNIT NAME. DEFAULT TO CANCEL. 
laOll 


Two plots are then displayed: one for the response of y sp i to a step in y spd i, and the other for y sp 2 to a 
step in y spd2 . 




STEP RESPONSES FOR NON-ZERO SET-POINT REG. 

input: set point command yspd < 2 ) 

AMPLITUDE = 1.00 

DEC 7, 1982 



2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 22.5 25.0 


With this low weighting on x\ 9 the system responses are obviously not acceptable, a fact also 
discernible from an examination of the closed-loop eigenvalues displayed by function 801. Thus the 
user requests another design iteration. 


function o 

TO COMPUTE FURTHER * ENTER NEXT FUNCTION N0S.<I3)p ONE PER LINE 


TO TERMINATE ENTER 999CLAST ENTRY MUST BE A RETURN) 
:203 
801 
301 
819 
604 


This time, the weight on xi is increased to 1.0, and the LQR and feedforward gains are 
recomputed. 


FUNCTION 203 

DISPLAY CONPAR BEFORE MAKING CHANGES? 
n 

ENTER CONTROL WTS QCp NN AND/OR PCINV (NAMELIST CONPAR) 
Seonpsr gc<1p1) = 10. Sem# 

Seonpar qc<1p1) - 10.0 Send 
£ CONPAR 

QC = 10. Op 2499*0.0 
NN = 250*0.0 

PC I NV = I.Op 5*0. Op 1.0- 18*0.0 
SEND 

ARE THERE ANYMORE CHANGESp Y OR N? 
n 


FUNCTION 801 

DESIGN A LINEAR QUADRATIC REGULATOR 


EIGENVALUES OF REGULATOR OR FILTER IN FN/ZETA FORM 
NAT FREQ (HZ) ZETA 

0.1392 1.000 

0.3027 0.5564 

SS = 


1 


3 


1 

2 

3 


7.206 

2.783 

0.9032 


2.783 0.9032 

1.950 0.6620 

0.6620 0.9392 


KC = 


1 2 3 


1 0.9032 0.6620 0.9392 

2 2.783 1.950 0.6620 


FUNCTION 301 
FORM A-BKC MATRIX 


FUNCTION 819 

FORM FEED FORWARD MATRIX FOR NON-ZERO SET POINT CONTROL 
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KFF 


1 


1 1,069 0,9412 

2 2.978 -0.3380 


FUNCTION 604 

OBTAIN AND PLOT SELECTED STEP RESPONSES 
FOR THE NON-ZERO SET POINT LINEAR REGULATOR 

STATE TRANSITION MATRIX OF LINEAR REGULATOR FOR TIME STEP 0.2400 


1 


2 


3 


1 0.9089 0.1832 

2 -0.5167 0.5573 

3 -0.8230D-01 -0.2934 


0.7571 D-02 
0 . 5555P-01 
0.7855 


FORCED RESPONSE MATRIX OF LINEAR REGULATOR FOR TIME STEP 0,2400 


1 


2 


3 


1 0,2315 0.2420P-01 

2 -0.6794D-01 0.1856 

3 -0.1475D-01 -0.3928D-01 


0. 6469 D -03 
0 . 7636D-02 
0.2138 


Now, the damping ratio of the complex eigenvalue pair has increased to 0.388, which is still too 
small to give acceptably damped responses. This fact is confirmed by the set -point responses, which 
are displayed next. 
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Actually, the 10-percent overshoot criterion is met, but y spl does not quite reach steady state in less 
than 5 seconds. Thus the user begins a third design iteration, choosing the weighting on xj to be 10 
this time. 


FUNCTION 0 

TO COMPUTE FURTHER , ENTER NEXT FUNCTION NOS* < 13), ONE PER LINE 
TO TERMINATE ENTER 999(LAST ENTRY MUST BE A RETURN) 

203 

801 

301 

819 

604 


FUNCTION 203 

DISPLAY CONPAR BEFORE MAKING CHANGES? 
n 

ENTER CONTROL WTS GC , NN AND/OR PCINV (NAMELIST CONPAR) 
SCONPAR 

QC= 1*0, 2499*0 » 0 
NN= 250*0.0 

PCI NV = 1*0, 5*0.0, 1.0, 18*0.0 
SEND 

ARE THERE ANYMORE CHANGES, Y OR N? 
n 


FUNCTION 801 

DESIGN A LINEAR QUADRATIC REGULATOR 


EIGENVALUES OF REGULATOR OR FILTER IN FN/ZETA FORM 
NAT FREQ (HZ) ZETA 

0.9524D-01 1.000 

0.2063 0.3879 


The preceding eigenvalues indicate that response time should now be acceptable as the slowest 
eigenvalue (Xi =0.1392 Hz = 0.875 rad/sec) should give a settling time (~4/Xi) of about 4.5 seconds. 


ss = 


i 


3 


1 1.285 0.6844 

2 0.6844 0.7483 

3 0.5240 0.3449 


0.5240 

0.344V 

0.7536 


KC = 


1 2 3 


1 0.5240 0.3449 0.7536 

2 0.6844 0.7483 0.3449 


FUNCTION 301 
FORM A-BKC MATRIX 


FUNCTION 819 

FORM FEED FORWARD MATRIX FOR NON-ZERO SET POINT CONTROL 
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KFF 


1 


1 0.6585 0 ♦ 7556 

2 0.7592 -0.6551 


FUNCTION 604 

OBTAIN AND PLOT SELECTED STEP RESPONSES 
FOR THE NON-ZERO SET POINT LINEAR REGULATOR 

STATE TRANSITION MATRIX OF LINEAR REGULATOR FOR TIME STEP 0.2400 


1 


3 


1 0.9574 0.2138 

2 -0.1549 0.7967 

3 -0 ♦ 8876D-0 1 -0.2785 


0*1 648D-0 1 

0.1292 

0.8125 


FORCED RESPONSE MATRIX OF LINEAR REGULATOR FOR TIME SI EP 0.2400 


1 


2 


3 


1 0.2356 0. 2673 D -01 

2 -0 ♦ 1901D-01 0.2164 

3 -0.1208D-01 -0 . 3521D-01 


0.1366 D- 02 
0 * 1 662D-0 1 
0.2177 


The step responses confirm that both set-point output transients are acceptable. 
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As a final check, the control responses are examined for this design to see that required control 
excursions are within limits (±5.0). To do this, the user requests function 101, which allows changes 
to be made in reference values. In particular, the user changes elements of transient response 
selection matrices MSYSP and MSPU so as to request both the control responses and the output 
responses, thus displaying the interaction between output 1 and input 2 and vice-versa. 


FUNCTION 0 

TO COMPUTE FURTHER, ENTER NEXT FUNCTION N0S.<I3>r ONE PER LINE 
TO TERMINATE ENTER 999CLAST ENTRY MUST BE A RETURN) 

101 

604 


FUNCTION 101 
CHANGES TO REFS 

DISPLAY REFS BEFORE MAKING CHANGES? 
n 

ENTER CHANGES 10 NAMELIST REFS (TSFTRr DT r Fir DELFr ZERMAX? AMPSPr AMF'SRr 
AMPICX, IFr ISPACEr IOUTr I ME AS r JINCr JINDr ITRMXr NCURUr L INLOG r MSPYr 
MSPYSPr MSPU r MSROLYr MSROLXr MICCLYr MICCLXr MICCLUr MICOLYr MICOLX) 


The program then displays the reference values and plots the requested transients. 


8 REFS mspysP<lrl>=0rlrllkSpysp(lr2) = lr0rinsPU<lrl> a lrlrmSFU<lf2> s: lrl Send 
8REFS 

TSFTR= 1.0 
DT= 0.240 
F I = 0 ♦ 1 0 D - 0 1 
DELF= 0.10D-01 
ZERMAX= 100.0 
AMPSP= 5* l. . 0 
AMPSR= 5*1.0 
A M P I C X = 50*1.0 
IF— 96 
ISPACE= J 
I0UT= 1 
I ME AS= 1 
J I NC= 1 
J I Nil** 1. 

I T R M X - 100 
NCURU- 2 
L I NLOG- 2 
MSPY- 250*0 

MSPYSPr Or lr 2*0, 1, 19*0 
MSPU= 2*1- 3*0 * 2*1. 18*0 
MSR0LY= 250*0 

MSR0LX= 2*1 r S *0 r 2*1 r 238*0 

M I CCL Y= lr 2499*0 

M I CCLX= 2500*0 

M I C C L U = 250*0 

MICOLY* lr 2499*0 

MI C0LX= lr 50* 0 r 1. 2448*0 

8 END 

ARE THERE ANYMORE CHANGES r Y OR N? 
n 


FUNCTION 604 

OBTAIN AND PLOT SELECTED STEP RESPONSES 
FOR THE NON-ZERO SET POINT LINEAR REGULATOR 

STATE TRANSITION MATRIX OF LINEAR REGULATOR FOR TIME STEP 0.2400 


1 


2 


3 


1 0.9089 0.1832 

2 -0.5167 0.5573 

3 -0 . 8230D-01 -0.2934 


0.7571D-02 
0 • 5555D-01 
0.7855 
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STEP RESPONSES FOR NON-ZERO SET-POINT REG. 

input: set point command tspd < 2 ) 

AMPLITUDE = 1.00 

DEC 7* 1982 



Note that there is some interaction between input 1 and output 2 but little between input 2 and 
output 1. 

Finally, the four control responses are displayed. These responses confirm that design criterion 3 
has been met. That is, both control signal absolute magnitudes remain less than 5 during a step 
response. The maximum excursion for U! is + 1.07 and that for u 2 is +2.98, both well below the 
design limits. 







After the responses are displayed, the user is again prompted for more requests, at which time the 
user terminates the program. The terminal session ends with the user requesting a printout of the 
output dataset (OUTxxx) that was generated during the present run for further off-line analysis if 
desired. 

FUNCTION 0 

TO CONFUTE FURTHER » ENTER NEXT FUNCTION NOS. (13) > ONE PER LINE 
TO TERMINATE ENTER 999CLAST ENTRY MUST BE A RETURN) 

999 

STORE THE N1 FOR THIS RUN? 
n 

CHCRW400 terminated: stop return 

Erint out33fprtsp=edit 

CZABD050 PRINT BSN=0498r 800 LINES 

off 

BOO? LOGOFF AT 14M0 ON 07/07/81 - CPU T I ME= 0.07 MINUTES. 

LOGICAL DISCONNECT » LOGON OR HANG UP 



Appendix D 

Terminal Output Options and Main PROCDEF 


This appendix includes (1) a list of those items included in the “standard” terminal output, (2) a 
list of those items included in the “extended” terminal output, and (3) the PROCDEF that is used to 
set up the program before it is run. 

Standard Terminal Output 

The following data will be displayed in the user’s terminal if the user has requested the functions 
that generate these data: 

(1) NAMELIST N1 

(2) NAMELIST REFS 

(3) NAMELIST CONPAR 

(4) NAMELIST ESTPAR 

(5) NAMELIST MATDAT 

(6) Open- and closed-loop eigenvalues 

(7) Kalman filter eigenvalues 

(8) Transfer function numerator and denominator polynomial coefficients 

(9) Transfer function gains 

(10) Maximum and average symmetry error for the SS and PP matrices 

(11) Positive-definiteness checks for the SS and PP matrices 

(12) Maximum element of and trace of the residual error matrix for the control and Kalman filter 
Riccati equations 

(13) Lyapunov error check data consisting of 

(a) Trace of residual 

(b) Normalized diagonal elements of error matrix 

(c) Trace of error 

(d) Trace of covariance 

(e) Ratio of trace of error to trace of covariance 

(14) Normalizing factors 

(15) Error messages 

Extended Terminal Output 

The extended terminal output consists of all of the standard terminal output plus 

(1) Input matrices 

(2) All eigenvectors and mode shapes 

(3) The ATOT, CTOT, DTOT, KCTOT, and HTOT matrices 

(4) The SS matrix 

(5) The PP matrix 

(6) The KC matrix 

(7) The KE matrix 

(8) The KFF matrix 

(9) All covariance matrices 

(10) Controllability check matrix 

(11) Observability check matrix (through H) 

(12) Observability check matrix (through C) 

(13) Transfer function gains and zeroes 

(14) All state-transition matrices 

(15) All forced-response matrices 

(16) Symmetry error matrix for the control Riccati equation 

(17) Residual error matrix for the control Riccati equation 

(18) Symmetry error matrix for the Kalman filter Riccati equation 



(19) Residual error matrix for the Kalman filter Riccati equation 

(20) Normalized system matrices: A, B, C, H, QQ, RRINV, D, DOUT, and CSP 

PROCDEF AESRUN 

The user invokes the following PROCDEF, AESRUN, before running the AESOP program. The 
purpose of the PROCDEF is to datadef all necessary libraries and required datasets. The PROCDEF 
requires one parameter (up to three characters), which is used to label all datasets that may be 
generated during the subsequent run of the AESOP program. 


AESRUN 

AESRUN 

AESRUN 

AESRUN 

AESRUN 

AESRUN 

AESRUN 

AESRUN 

0000000 

0000100 

0000200 

0000300 

0000400 

0000500 

0000600 

0000700 

AESRUN 

0000800 

AESRUN 

0000900 

AESRUN 

0001000 

AESRUN 

0001100 

AESRUN 

AESRUN 

AESRUN 

0001200 

0001300 

0001400 

AESRUN 

0001500 

AESRUN 

0001600 

AESRUN 

0001700 


PROCDEF AESRUN 
PARAM $1 
ERASE OUTS1 

DDEF FT06F001,VS,OUT$1,DCB = (RECFM = V,LRECL = 132), RET = T 
DISPLAY 'OUT$l IS DDEFD TO THE HI-SPEED PRINTER' 

DDEF FT08F001 ,. VS,CG$1 ; DISPLAY 'CG$1 IS DDEFD TO IO UNIT 8' 
DDEF FT09F00 1 , VS ,EG$ 1 ; DISPLAY 'EG$1 IS DDEFD TO IO UNIT 9' 
DDEF FT10F001,VS,PFRUZ$1; DISPLAY 'PFRUZSl IS DDEFD TO IO 
UNIT 10' 

DDEF FT11F001,VS,PFRUY$1; DISPLAY 'PFRUY$1 IS DDEFD TO IO 
UNIT 1 1 ' 

DDEF FT12F001 , VS,PFRWZ$1 ; DISPLAY 'PFRWZ$1 IS DDEFD TO IO 
UNIT 12' 

DDEF FT13F001 , VS,PFRWY$1 ; DISPLAY 'PFRWYSl IS DDEFD TO IO 
UNIT 13' 

DDEF FT14F001 , VS,CFR$1 ; DISPLAY 'CFRS1 IS DDEFD TO IO UNIT 
14' 

DDEF FT15F001 , VS,PP$1 ; DISPLAY 'PP$1 IS DDEFD TO IO UNIT 15' 
DDEF FT16F001,VS,SS$1; DISPLAY 'SS$1 IS DDEFD TO IO UNIT 16' 
DDEF FT17F001 , VS,FFG$1 ; DISPLAY 'FFGSl IS DDEFD TO IO UNIT 
17' 

DDEF RUN1, VP .GRAPHICS, OPTION = JOBLIB; DISPLAY 'RUN1 IS 
LIBRARY GRAPHICS' 

DDEF RUN2, VP, AESLIB, OPTION = JOBLIB: DISPLAY 'RUN2 IS 
LIBRARY AESLIB' 

LOAD BLOCKA9; LOAD AESOP 


Eleven VS datasets that are to contain program input or output are datadeffed by this PROCDEF. 
Table I lists these datasets plus four others that might be datadeffed by the user during the course of 
running the program. 
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Appendix E 
Flow Chart 


This appendix contains a subroutine flow chart of the numbers are listed (in parentheses) below the name of the 
AESOP program in “tree” form. Where specific func- subroutine, 
tions are performed in a subroutine, the function 


AESOP 



— MATPRT 


—MATPRT 


— EIGEN — 

( 401 ) 

—ARRAY 
— HSBG 
— EIGQR 
— SC ALEA 
—CONDI — 


— EGVCTR — 

( 402 ) 

— MATPRT 
—ARRAY 
— FACTR 
— PRMUTE 

— M0D5HP — 

( 402 ) 

—MATPRT 


— 5CALEA 
— REDU 
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— CTBL 
(403) 


— AES500- 
(503,506, 
(509,512, 
(525) 


--AES6 0 0- 


— AES7 00- 


— OBSBL 
(403) 


--RESI 

(403) 


— NRML 

(404) 

— UNRML 

(405) 


— MATPRT 
— MXINV 

— MATPRT 

—MATPRT 


—MATPRT 


) 

) 


— PREREQ 

— PRSPN5 — 
(501 , 504, ) 
(507,510, ) 
(513,515, ) 
(517,519, ) 
(521,523) 


— FRQP — 

— BOLLIN— 


— FRPOL Y 

— DA V I SO 
— DAH5KY-- 


— BODE —| 
(502,505, ) | 
(508,511, ) | 
(514,516, ) 
(518,520, ) 
(522,524) | 


-PLOTSUBS 


— PREREQ 

—MATPRT 

— DSCRT 
(601,602,) 

(603.604) 

— STP 

(601.604) 


— ICR5P — 
(602,603) 


| -PLOTSUBS 


-PLOTSUBS 


— PREREQ 


—MATPRT 


— ZEROES— 
(701,702, ) 
(703,704, ) 
(705) 


—ARRAY 


— HSBG 


— EIGQR 
— CONDI — 


— 5CALEA 


— POLMPY 


— REDU 



— GAIN 
(701,702, ) 
(703,704, ) 
(705) 


— AES80 0 — 
(802 ,806 , ) 
(808,810, ) 
(814,816) 


— PREREQ 

— MATPRT 

— CONTRl — 
(801) 


—MATPRT 
— RICSS - 


—MATPRT 

— ARRAY 
— MXINV 

— ORDER 
— EIGQR 
— 5CALEA 

— HSBG 

— EGCK —| 


—MODSHP — 


—MATPRT 


--MATPRT 
— ARRAY 
— FACTR 
— PRMUTE 


— EIGEN — 
(803,805, ) 
(811,812, ) 
(813) 


—ARRAY 


—HSBG 
— EIGQR 
— SCALEA 
— CONDI — 


—SCALEA 

--REDU 


— EGVCTR — 
( 804 ) 


— MATPRT 
—ARRAY 
— FACTR 
—PRMUTE 


—MODSHP — 
(804) 


—MATPRT 


— RICCHK— 
(807,815) 


—MATPRT 


— ESTMAT— 
(809) 


—MATPRT 
— RICSS — 


—MATPRT 
— ARRAY 


—MXINV 



— COVAR 
(817) 


— LYPCK 
(818) 


— MXINV 
(819) 


—ORDER 
— EIGQR 
— SCALEA 
— HSBG 
— EGCK 
— MODSHP— | 

j — MATPRT 
— EGVCTR— j 

--MATPRT 
— ARRAY 
j — FACTR 
| — PRMUTE 


— MATPRT 
— ARRAY 
— LAPNV — j 

1--MXADD 

I — MXMIT 
J — MXTRA 
—MXINV 
| — DSCA 


— MATPRT 
— ARRAY 
--LAPNV — 


— MXADD 
— MXMLT 
— MXTRA 
—MXINV 
— DSCA 


— AE5900-- 
(901,902, ) 
(903,909) 


— UZR90 1 
— UZR902 
— U2R903 
— UZR909 



Appendix F 
Prerequisite Table 


This appendix contains the information used by the 
AESOP program to make prerequisite checks. These 
checks are performed before each AESOP function is to 
be executed to see that the necessary function or com- 
bination of functions has been performed prior to 
execution of the present function. Table XII lists this 
prerequisite check information in the following form: 
The specific prerequisite function or combinations of 


functions are listed across the top, and each function 
whose prerequisites are to be checked is listed on the left. 
A checkmark appears in a row to indicate a prerequisite. 
Multiple checkmarks in a row indicate that the corre- 
sponding prerequisites are to be logically “ANDed.” For 
example, the logical prerequisite statement for function 
303 is (201 OR 202) AND (205 OR 801) AND (206 OR 
809). 
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TABLE XII. - PREREQUISITES FOR AESOP FUNCTIONS 





TABLE XII 


Concluded 


Function Prerequisite function 

None 201 205 206 207 208 301 302 303 205, 401 ' 402 404 501 504 507 510 513 ' 515 517 519 521 523 801 803 809 817 819 

or or or or or 206, or 

202 801 809 801 809 801, 2U9 

or 
809 


522 

523 

524 

525 


X 

X 


601 X 

602 X 

603 

604 X 


X 

X 


701 X 

702 X 

703 X 

704 X 

705 


X 


801 

802 

803 

804 

805 

806 

807 

808 

809 

810 
811 
812 

813 

814 

815 

816 

817 

818 

819 

820 

901 

902 

903 

904 


X 


X 


X 


X 

X 

X 

X 


X 


X 

X 


X 

X 


X 


X 


X 

X 


X 

X 


X 
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